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1 Introduction 

The basic idea of density functional theory (DFT) is to describe an interacting 
many-particle system exclusively and completely in terms of its density. The for- 
malism rests on two basic theorems: 
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I. Every observable quantity can be calculated, at least in principle, from the 
density alone, i. e. each quantum mechanical observable can be written as a 
functional of the density. 

II. The density of the interacting system of interest can be obtained as the density 
of an auxiliary system of non-interacting particles moving in an effective local 
single-particle potential, the so-called Kohn Sham potential. 

In the original work of Hohcnbcrg and Kohn (HK) [1] and Kohn and Sham (KS) 
[2] these theorems were proven for the ground-state density of static many-body 
systems. On the basis of these theorems, DFT has provided an extremely success- 
ful description of ground-state properties of atoms, molecules and solids [3, 4, 5]. 
The quality of approximations for the Kohn-Sham potential has steadily improved 
over the years and the currently best functionals yield ground-state properties in 
very close agreement with configuration interaction results [6]. Excited-state prop- 
erties, however, are notoriously difficult to calculate within the traditional density 
functional framework and time-dependent phenomena are not accessible at all. 

Time-dependent density functional theory (TDDFT) as a complete formalism 
[7] is a more recent development, although the historical roots date back to the 
time-dependent Thomas-Fermi model proposed by Bloch [8] as early as 1933. The 
hrst and rather successful steps towards a time-dependent Kohn-Sham (TDKS) 
scheme were taken by Peuckert [9] and by Zangwill and Soven [10]. These authors 
treated the linear density response of rare-gas atoms to a time-dependent external 
potential as the response of non-interacting electrons to an effective time-dependent 
potential. In analogy to stationary KS theory, this effective potential was assumed 
to contain an exchange-correlation (xc) part, v xc (r,t), in addition to the time- 
dependent external and Hartree terms: 

v a (r,t)=v(r,t) + Jd 3 r' ^^+v xc (r,t) . (1) 

Peuckert suggested an iterative scheme for the calculation of v xc , while Zangwill 
and Soven adopted the functional form of the static exchange-correlation potential 
in local density approximation. Significant steps towards a rigorous foundation of 
time-dependent density functional theory were taken by Deb and Ghosh [11]- [14] 
and by Bartolotti [15] -[18] who formulated and explored HK and KS type theorems 
for the time-dependent density. Each of these derivations, however, was restricted 
to a rather narrow set of allowable time-dependent potentials (to potentials periodic 
in time in the theorems of Deb and Ghosh, and to adiabatic processes in the work 
of Bartolotti) . A general proof of statements I and II above for the time-dependent 
density was given by Runge and Gross [7] . A novel feature of this formalism, not 
present in ground-state density functional theory, is the dependence of the respective 
density functionals on the initial (many-particle) state. A detailed description of 
the time-dependent density functional formalism will be presented in section 2. 
The central result is a set of TDKS equations which are structurally similar to the 
time-dependent Hartree equations but include (in principle exactly) all many-body 
effects through a local time-dependent exchange-correlation potential. In section 2 
we focus on the motion of electrons only. In many experimental situations, however, 
the nuclear motion and its coupling to the electronic motion is important as well. If, 
for example, a molecule is placed in the focus of a strong laser, the electric field of 
the laser can either couple directly to the nuclei (in the infrared frequency regime) or 
the coupling to the electrons can lead to photoionization with subsequent Coulomb 
explosion (dissociation) of the molecule. To deal with situations of this type a DF 
formalism for the coupled system of electrons and nuclei is developed in section 
3. The central result is a set of coupled TDKS equations for the electrons and for 
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each nuclear species. While sections 2 and 3 exclusively deal with time-dependent 
electric fields, magnetic effects will be considered in section 4. Both the ordinary 
Zccman coupling and the coupling of magnetic fields to the orbital currents will be 
included. 

To date, most applications of TDDFT fall in the regime of linear response. The 
linear response limit of time-dependent density functional theory will be discussed 
in section 5.1. After that, in section 5.2, we shall describe the density-functional 
calculation of higher orders of the density response. For practical applications, 
approximations of the time-dependent xc potential are needed. In section 6 we 
shall describe in detail the construction of such approximate functionals. Some 
exact constraints, which serve as guidelines in the construction, will also be derived 
in this section. Finally, in sections 7 and 8, we shall discuss applications of TDDFT 
within and beyond the perturbative regime. Apart from linear response calculations 
of the photoabsorbtion spectrum (section 7.1) which, by now, is a mature and 
widely applied subject, we also describe some very recent developments such as the 
density functional calculation of excitation energies (section 7.2), van der Waals 
forces (section 7.3) and atoms in superintense laser pulses (section 8). 

2 Basic formalism for electrons in time-dependent 
electric fields 

2.1 One-to-one mapping between time-dependent potentials 
and time-dependent densities 

Density functional theory is based on the existence of an exact mapping between 
densities and external potentials. In the ground-state formalism [1], the existence 
proof relies on the Rayleigh-Ritz minimum principle for the energy. Straightforward 
extension to the time-dependent domain is not possible since a minimum principle 
is not available in this case. The proof given by Runge and Gross [7] for time- 
dependent systems is based directly on the Schrodinger equation (atomic units are 
used throughout): 

i^(t) = H(tmt)- (2) 

We shall investigate the densities n(r, t) of electronic systems evolving from a fixed 
initial (many-particle) state 

*(to) = *o (3) 
under the influence of different external potentials of the form 

JV 

V(t) = 5>( ri ,t). (4) 

i=l 

In the following discussion, the initial time to is assumed to be finite and the po- 
tentials are required to be expandable in a Taylor scries about t . No further 
assumptions concerning the size of the radius of convergence are made. It is suf- 
ficient that the radius of convergence is greater than zero. The initial state is 
not required to be the ground state or some other eigenstate of the initial potential 
v(r,to) = vq(t). This means that the case of sudden switching is included in the 
formalism. On the other hand, potentials that are switched-on adiabatically from 
t = — oo are excluded by the Taylor-expandability condition because adiabatic 
switching involves an essential singularity at to = — oo. 

Besides an external potential of the form (4) , the Hamiltonian in Eq. (2) contains 
the kinetic energy of the electrons and their mutual Coulomb repulsion: 

H(t) = f + U + V(t) (5) 
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with 

T = Y.-\^i (6) 

i=l z 



and 

AT 



2^|r 4 -r J | 



(7) 



With these preliminaries, we can formulate the following Hohenbcrg-Kohn-type 
theorem: The densities n(r, t) and n'(r, t) evolving from a common initial state 
= ^(to) under the influence of two potentials v(r,t) and v'(r,t) (both Taylor 
expandable about the initial time to) are always different provided that the poten- 
tials differ by more than a purely time-dependent (r-independent) function: 1 

v(r,t)^v'(r,t)+c(t). (8) 

To prove this theorem, we use the condition that the potentials v and v' can be 
expanded in Taylor series: 

oo 

v(r,t) = J2vV k m-t ) k , (9) 

oc 

v'(r,t) = ^^ k (v)(t-t Q ) k . (10) 

k=0 

Equation (8) is equivalent to the statement that for the expansion coefficients v k (r) 
and v' k (r) there exists a smallest integer k > such that 



Qk 

v k (r)-v' k (r) = w (v(r,t)-v'(r,t)) 



7^ const. (11) 

t=t 



From this inequality we prove in a first step that the current densities 

j(r,t) = <*(t)U P (r)|tt(t)> (12) 

and 

j'(r,t) = <*'(t)|j p (r)|*'(t)) (13) 
are different for different potentials v and v'. Here, 

1 N 

JpW = 2i £ (V'A T - T J) + 5 ( T - r i) V *,) ( 14 ) 

i=i 

is the usual paramagnetic current density operator. In a second step we shall show 
that the densities n and n' are different. 

Using the quantum mechanical equation of motion for the expectation value of 
an operator Q(t), 

!<*(t)|Q(t)|*(t)) - <*(t)| ( ^ - i[Q(t),#(*)] ) !*(*)> , (15) 



1 If u and i/ differ by a purely time-dependent function, the resulting wave functions \P(t) and 
\l/'(t) differ by a purely time-dependent phase factor and, consequently, the resulting densities p 
and p' are identical. This trivial case is excluded by the condition (8), in analogy to the ground- 
state formalism where the potentials are required to differ by more than a constant. 
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we obtain for the current densities: 



^j(r,i) = |<*(t)|jp(r)|*(t)) = -i<*(t)|[j p (r), J ff(t)]|*(t)) 

|j'(r,i) = |<*'(t)|jp(r)|*'(t)) = -i<*'(t)|[j p (r),ff'(t)]|*'(t)) 
Since * and evolve from the same initial state 

tf (to) = *'(*o) = *o 

we can write 



(16) 
(17) 

(18) 



Si 



(j(r,t)-j'(r,t)) 



t=t 



= -i(* \\j p (r),H(t )-H'(toWo) 
= -n (r)V{v{r,t )-v'{r,t )) 



with the initial density 



n (r) = n(r,i ) . 



(19) 



(20) 



If the condition (11) is satisfied for k — the right-hand side of (19) cannot vanish 
identically and j and j' will become different infinitesimally later than to- If (H) 
holds for some finite k > we use Eq. (15) (k + 1) times and obtain after some 
algebra: 



fe+i 



(j(r,t)-j'(r,t)) 



-n (r)Vw fc (r) ^ 



t=t 



with 



w k (r) = 

Once again, we conclude that 



(v{r,t)-v'(r,t)) 



t=t 



j(r,i)#j'(r,i) 



(21) 
(22) 

(23) 



provided that (11) holds for v and v'. To prove the corresponding statement for the 
densities we use the continuity equation 



d_ 

di 



(n(r,t)-n'(r,t)) = -V-(j(r,t)-j'(r,t)) 



and calculate the (k + l)th time derivative of Eq. (24) at t = t a : 

k+2 



(n(r,t)-n'(rt)) 



V • (no(r)Vw fe (r)) . 



(24) 



(25) 



t=t 



In order to prove that the densities n(r, t) and n'(r, t) will become different infinites- 
imally later than to, we have to demonstrate that the right-hand side of Eq. (25) 
cannot vanish identically. This is done by reductio ad absurdum: Assume 

V-(n (r)V Wfe (r))=0 (26) 

and evaluate the integral 
d 3 rn (r)[Vw fe (r)] 2 

= - Jd 3 rw k (r)V • (n (r)V«7 fe (r)) + £ dS ■ (n (r)w k (r)Vw k (r)) , (27) 
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where we have used Green's theorem. The first integral on the right-hand side 
of (27) vanishes by assumption. For physically realistic potentials (i. e. potentials 
arising from normalizable external charge densities) , the surface integral vanishes as 
well, because for such potentials the quantities u>fc(r) fall off at least as 1/r. Since 
the integrand on the left-hand side is non-negative one concludes that 

n (r)[V^(r)] 2 = (28) 

in contradiction to u>fc(r) 7^ const. This completes the proof of the theorem. We 
mention that more general potentials may also be considered. The precise conditions 
have been formulated in [19]. 

We note in passing that the right-hand side of Eq. (25) is linear in ■ Conse- 
quently, the difference between n(r,i) and n'(r,t) is non-vanishing already in first 
order of v(r,t) — v'(r,i). This result will be of importance in section 5 because it 
ensures the invertibility of linear response operators. 

By virtue of the 1-1 correspondence established above (for a given ^o), the 
time-dependent density determines the external potential uniquely up to within 
an additive purely time-dependent function. The potential, on the other hand, 
determines the time-dependent wave function, which can therefore be considered 
as a functional of the time-dependent density, unique up to within a purely time- 
dependent phase a(t): 

*(t) = e - M (*)*[n](t) . (29) 

As a consequence, the expectation value of any quantum mechanical operator Q(t) 
is a unique functional of the density: 

Q[n](t) = <*[n](t)|$(t)|*[n](t)> . (30) 

The ambiguity in the phase cancels out. As a particular example, the right-hand side 
of Eq. (16) can be considered as a density functional which depends parametrically 
on r and t: 

P[n](r,t) = -i<*[n](t)|[j p (r),H(t)]|*[n](t)) . (31) 

This implies that the time-dependent particle and current densities can always 
be calculated (in principle exactly) from the following set of "hydrodynamical" 
equations: 

^n(r,i) = -V-j(r,f) (32) 

^j(r,i)=P[n](r,t) . (33) 
In practice, of course, the functional P[n] is only known approximately. 



2.2 Stationary-action principle 

The solution of the time-dependent Schrodinger equation (2) with initial condition 
(3) corresponds to a stationary point of the quantum mechanical action integral 

A = j\t^(t)\i^-H(t)Mt))- (34) 

Since there is a 1-1 mapping between time-dependent wave functions, ^(t), and 
time-dependent densities, n(r, t), the corresponding density functional 

A = j\t <*[n](*)K! - H(t)Mn](t)) (35) 
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must have a stationary point at the correct time-dependent density (corresponding 
to the Hamiltonian H(t) and the initial state ^o)- Thus the correct density can be 
obtained by solving the Euler equation 

SA[n] -0 (36) 



dn(r, t) 

with appropriate boundary conditions. The functional A[n] can be written as 

A[n] = B[n] - jf dt J d 3 r n(r, t)v(r, t) (37) 

with a universal (^-dependent) functional B[n], formally defined as 

B[n] = j\t (*[n](t)|i| - t - U\*[n]{t)) . (38) 

On the exact level, the hydrodynamical equations (32, 33) and the variational equa- 
tion (36) are of course equivalent. The functionals P[n], A[n], B[n] are well-defined 
only for v-representable densities, i. e. for densities that come from some time- 
dependent potential satisfying Eq. (9). In view of this, a Levy-Lieb-type [20, 21, 22] 
extension of the respective functionals to arbitrary (non-negative, normalizable) 
functions n(r, t) appears desirable. Two different proposals of this type have been 
put forward so far [23, 24]. 



2.3 Time-dependent Kohn-Sham scheme 

The 1-1 correspondence between time-dependent densities and time-dependent po- 
tentials can be established for any given interaction U , in particular also for U = 0, 
i. e. for non-interacting particles. Therefore the external potential v s [n](r,t) of a 
non-interacting system reproducing a given density n(r, t) is uniquely determined. 
However, the 1-1 correspondence only ensures the uniqueness of v s [n] for all v- 
representable densities but not its existence for an arbitrary n(r,t). In order to 
derive a time-dependent KS scheme we have to assume, similar to the static case, 
non- interacting w-representability, i.e., we have to assume that a potential v s exists 
that reproduces the time-dependent density of the interacting system of interest. 
Under this assumption, the density of the interacting system can be obtained from 

N 

n(r,t) = ^|^-(r,t)| 2 (39) 

with orbitals ipj (r, t) satisfying the time-dependent KS equation 

d /V 2 \ 

i-Q t <Pj(r,t) = \-—+v.[n](T,t))<pj(T,t) . (40) 

Usually, the single-particle potential v s is written as 

v s [n](r,t)=v{r,t) + Jd 3 r'^^-+v xc [n]{r,t) , (41) 

where v(r, t) is the external time-dependent field. Equation (41) defines the time- 
dependent xc potential. In practice, this quantity has to be approximated. As in the 
static case, the great advantage of the time-dependent KS scheme lies in its compu- 
tational simplicity compared to other methods such as time-dependent Hartree-Fock 
or time-dependent configuration interaction [25]- [32]. One has to emphasize that, 
in contrast to time-dependent Hartree-Fock, the effective single-particle potential 
v s is a local potential, i.e., a multiplicative operator in configuration space. 
A few remarks are in order at this point: 
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(i) An important difference between the ordinary ground state density functional 
theory and the time-dependent formalism developed above is that in the time- 
dependent case the 1-1 correspondence between potentials and densities can 
be established only for a fixed initial many-body state Consequently, 
the functionals P[n], A[n] and B[n] implicitly depend on ^o- In the same 
way, v s [n] and w xc [n] implicitly depend on the initial KS Slater determinant. 
The formalism provides no guideline of how to choose the initial KS orbitals 
(fij(r,to) as long as they reproduce the initial interacting density no corre- 
sponding to *o- In general, there exist infinitely many Slater determinants 
reproducing a given density [33, 34]. From a formal point of view there is 
no problem with that; any choice of initial orbitals (pj(r,to) reproducing the 
initial interacting density no will do the job because the dependence of v s [n] 
on the initial state is such that the interacting density will be reproduced 
in each case. In practice, however, the dependence on the initial state is a 
nuisance. Of course one would prefer to have functionals of the density alone 
rather than functionals of n(r,t) and ^>o- One has to emphasize, however, 
that for a large class of systems, namely those where both and the initial 
KS Slater determinant are non-degenerate ground states, P[n] and v s [n] are 
indeed functionals of the density alone. This is because any non-degenerate 
ground state ^f is a unique functional of its density n (r) by virtue of the 
traditional HK theorem. In particular, the initial KS orbitals are uniquely 
determined as well in this case. 

(ii) We emphasize that the KS scheme does not follow from the variational prin- 
ciple. Incidentally, the same statement holds true in the static case as well. 
The KS scheme follows from the basic 1-1 mapping (applied to non-interacting 
particles) and the assumption of non-interacting ^-represent ability. The varia- 
tional principle yields an additional piece of information, namely the equation 

where A xc is the xc part of the action functional, formally defined by 

Ac N = B s [n] B[n] \ f\t J d\ J dV <^<^ 1 ) . (43) 

Here B s [n] is the non-interacting analogue of the functional B[n], i. e. , 

B B [n] = J\t <*[n](t)|i| - f\<S>[n]{t)) (44) 

where $[n](t) is the unique time-dependent Slater determinant corresponding 
to the density n. 

(iii) The current density 

1 N 

j(r,i) = -^(^(r,t)V^ fc (r,t)- Vfc (r,t)V^(r,t)) (45) 
1 fe=i 

following from the TDKS orbitals is identical with the true current density 
of the interacting system at hand. In order to prove this statement we recall 
that the first part of the Runge-Gross proof described above establishes a 1-1 
mapping, v(r, t) <-» j(r, t), between external potentials and current densities of 
interacting particles and likewise, for U = 0, a 1-1 mapping, v s (r, t) <-> j s (r, t), 
between external potentials and current densities of nonintcracting particles. 
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Making once again the assumption of non-interacting v-representability of the 
interacting current density j(r, t), one can establish an alternative "current- 
density version" of the KS scheme, 

= (-iv 2 +5 s [j](r,t)) ^.(r,i) (46) 

1 N 

j( r ' *) = 2i S ( ^ (r ' *) V ^( r - *) - <^( r ' *) V #E(r. *)) ( 47 ) 
fe=i 

whose solution reproduces the current density, j, of the interacting system of 
interest. To prove Eq. (45), we show that the solutions <pj(r,t) of (46,47) are 
in fact identical with the solutions ipj of the ordinary TDKS scheme (39-41). 
To this end we prove that the density 

N 

n(r,f)=^|^(r,i)| 2 (48) 

k=l 

is identical with the density resulting from (39-41): 

n(r, t) = n(r, t) . (49) 

Then the uniqueness of the potential v s [n] reproducing n(r, t) implies that 
v s (r,t) = v s (r,t) so that the solutions of (46) and (40) are identical. In order 
to prove Eq. (49) we observe that the full many-body Schrddinger Eq. (2) 
implies the continuity equation 

^ = -V.J(r,t) (50) 
while the Schrddinger equation (46) implies the continuity equation 

^ = -V.j(M). (51) 

Comparing (50) and (51) wc find that n(r,t) and n(r, t) can differ at most by 
a time-independent function r){r) so that, at the initial time to, 

n(r,to) = n(r,i )+fj(r). (52) 

Hence, if the initial orbitals are chosen to be identical, 

(fik(r,t ) = <fik(r,t ) k = l,...,N, (53) 

it follows that »j(r) = and Eq. (49) will be satisfied for all times. It remains 
to be shown that the choice (53) is always possible. This is not obvious 
a priori because, by construction, the orbitals <pfc( r >*o) must reproduce the 
initial density n(r, to) while the orbitals <£fc(r, t ) must yield the initial current 
density j(r, t )- In order to prove that the choice (53) is possible we show 
that a given density n (r) and a given current density jo(r) can always be 
simultaneously reproduced by a single Slater determinant 

$(r 1 ,...,r JV ) = -^=det{^(r < )}. (54) 

This can be shown with a current-density generalization [24] of the so-called 
Harriman construction [33]. Here we reproduce the construction for one spa- 
tial dimension. The three dimensional case can be treated in analogy to 
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Ref. [34]. Given the densities no (a;) and jo(x) we define the following func- 
tions 



q(x) :-- 



2tt 



x) := f 

J a 



dx' no(x') 
n (x') 



dx 



(55) 
(56) 



so that 



and 



dq(x) 2ir . 



(57) 
(58) 



ds(x) _ j (x) 
dx no(x) 

In terms of these quantities and the particle number N we define the single- 
particle orbitals 



4> k (x) := ] j^Me^ k ^+ s ^-^) 



(59) 



where k is an arbitrary integer while M is a fixed integer to be determined 
below. The functions : k integer} form a complete and orthonormal set 
(see e. g. [4]). Constructing the Slater determinant (54) from these orbitals it 
is readily verified that 

N 



n (x) 



(60) 



and 



.7 = 1 x ' 

ds n (x) dq 



n (x) 



dx 



N dx 




(61) 



Hence, by virtue of Eq. (58), Eq. (61) reproduces the given current density 
jo(x) if M is chosen equal to k ^j. This completes the proof. 



3 Motion of the nuclei 

3.1 Quantum mechanical treatment of nuclear motion 

The formalism developed so far is adequate whenever the motion of the atomic 
nuclei can be neglected. Then the electron-nucleus interaction only enters as a 
static contribution to the potential v(rt) in Eq. (41). This is a good approximation 
for atoms in strong laser fields above the infrared frequency regime. When the 
nuclei are allowed to move, the nuclear motion couples dynamically to the electronic 
motion and the situation becomes more complicated. 

In this section we shall describe a TDDFT for systems consisting of N electrons 
and Na nuclei of charge Za and mass Ma (in a.u.), A = 1, . . . , K. K is the number 
of different nuclear species. Let R-Aa be the configuration space vector of the ath 
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nucleus of species A. Then the complete system of electrons and nuclei is described 
by the Schrodinger equation 

d r - 

i — ^{r 1 ...r N ,{R Aa },t) = |# e (ri . . . rjv, t) + H n ({TL Aa }, t) 

+£r en (ri...rj V ,{R j4 a}) *(ri...rjv,{RAa},i) (62) 
with the electronic Hamiltonian, 



( 1 \ 1 1 



(63) 



the nuclear Hamiltonian 



K N A , , K N A K N B 

*» = E E f-^vL a + v&dw)) + JEEEE " 1 i; 



71 7 A=l a = l B=l 13=1 1 ^" Dpl 



(Aa)#(B/3) 

(64) 



and the electron-nucleus interaction 

N K N A 



^-EEE^n- (^) 

Based on an extension [35] of the Runge-Gross theorems described in sect. 2 to 
arbitrary multicomponent systems one can develop [36] a TDDFT for the coupled 
system of electrons and nuclei described above. In analogy to sections 2.1 - 2.3, 
one can establish three basic statements: First of all, there exists a rigorous 1-1 
mapping between the vector of external potentials and the vector of electronic and 
nuclear densities, 

(v mt (r,t);V£ ct (B.,t),...,V£ t (B.,t)) ^ (n(r,t);ni(R,t),...,n*(R,t)) • (66) 

Once again, this 1-1 correspondence is valid for a fixed initial many-body state 
^(ri, . . . rjv, {R^iq,}; to). Besides this HK-type statement, one can derive a stationary- 
action principle and a set of coupled TDKS equations for electrons and nuclei. The 
latter read as follows: 

i^Vj(r,t)= ^-^V2+w s [n,{n B }](r,^^(r,t), (67) 

j = l,...,N 

i^ Aa (K, t) = ^-—L V^, + V s A [n, {n B }](R, t)j V>Aa(R, t) , 

A = 1,...,K; a = l,...,N A 

with the nuclear densities 

N A 

n A (R,t) = J2n Aa (R,t), n Aa (R,t) = \ijj Aa {H,t)\ 2 , (69) 

a=l 

the electronic density 

N 

n(r,i)=E|^(r,i)| 2 , (70) 



(68) 
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and the KS potentials 



v a [n, {n B }](r,t) = v cxt (r, t) + d d r 



,3„/ n ( v ',t) 



-E fd'R 2 ^^ +v xc [n,{n B }](r,t), (71) 



n(r',t) 



V s A [n, {n B }](R, t) = V c A xt (Tl, t)-Z A J dV ^ 

+Z A J2 jd*R' ZB ^ B ^ ] + V A [n,{n B }](R,t). (72) 

The xc potentials in (71) and (72) are formally given by functional derivatives of 
the xc part of the quantum mechanical action functionals: 

-('■*)- '-^w 1 (73) 

In molecules and clusters, genuine exchange (as well as correlation) among iden- 
tical nuclei is very small because, at typical internuclear separations, the overlap 
of nuclear wave functions is rather small. 2 However, the exact xc functional also 
contains self-exchange contributions which are not small and which cancel the self- 
interaction terms contained in the Hartree potentials in Eqs. (71) and (72). Hence 
it will be a very good approximation to represent by the self-exchange terms 
alone. This leads to 

l3 j n(r',t) 



V s Aa [n,{n B }](R,t) = V e A t (R,t) - Z A J d 3 r' ^ 

B = l J ' ' (3 = 1 •> ' ' 



(75) 



Note that, within this approximation, the nuclear KS potential depends on the state 
ip Aa it acts on. This is analogous to the SIC scheme of Perdew and Zunger [37]. 

Clearly, a complete numerical solution of the coupled KS equations (67, 68) for 
electrons and nuclei will be rather involved. Usually only the valence electrons need 
to be treated dynamically. The core electrons can be taken into account approxi- 
mately by replacing the electron-nucleus interaction (65) by suitable pseudopoten- 
tials and by replacing the nuclear Coulomb potential in Eq. (64) by the appropriate 
ionic Coulomb potential [38]. This procedure reduces the number of electronic KS 
equations and hence the numerical effort considerably. 



3.2 Classical treatment of nuclear motion 

Further simplification is achieved by treating the nuclear motion classically. Numer- 
ical schemes of this type have been derived in various ways [38, 39, 40, 41, 42, 43]. 
In this section we shall use the multicomponent formalism developed in section 3.1 
as a starting point to derive classical equations of motion for the nuclei. 

Applying Ehrenfest's theorem to the nuclear KS equation (68), the classical 
trajectory 

R c ia SS W = (V>Aa(*)|R|Vu«(i)) = J d 3 R Kn Aa (R,t) (76) 

2 In atomic scattering processes this is not necessarily the case. 
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of the ath nucleus of species A satisfies the equation of motion 

M A ^^{t) = Y Aa {t) (77) 

where the force is given by 

F Aa (t) = -(*PA a (t)\V R V s A \tp Aa (t)) . (78) 

Since the TDKS equations (67) - (72) reproduce the exact nuclear densities, Eq. (76) 
yields the exact classical trajectory whenever species A contains only one nucleus. 
When species A contains more than one nucleus we have a system of indistinguish- 
able particles and then, strictly speaking, the trajectories of single nuclei cannot 
be told apart: Only the total density n A (R, t) and hence the center-of-mass tra- 
jectory of species A can be measured. In this case, trajectories of single nuclei can 
be defined by Eq. (76) within some effective single-particle theory TDKS theory is 
particularly suitable for this purpose since the TDKS partial densities n Aa lead to 
the exact total density n A ■ 

Employing the approximation (75) for the nuclear KS potential the force (78) 
on the nucleus (Act) simplifies to 

F Aa (t) = - Jd 3 Rn Aa (R,t) V R V e i t (R,t) - Z A J rf 3 rV R ^ 

B = l •> ' I 0=\ •> 



|R-R' 



(79) 



In many cases, the nuclear densities n Aa (R, t) will be rather narrow functions, with 
a strong peak at the classical trajectory R^ ss (i). In such a situation, integrals 
of the form fd 3 Rn Aa (TL, t)G(R) are well represented by Taylor-expanding G(R) 
around the classical trajectory. This leads to 



/ 



d 3 Rn Aa (K,t)G(K) = G(R c A T(t)) + 0(R - R^D • (80) 

The first-order term vanishes due to the definition (76) of the classical trajectory. 
Neglecting terms of second and higher order is equivalent to replacing the nuclear 
densities by (5-functions: 

n Aa (H,t) = <5(R-Rtf s (t)). (81) 
In this way the Newton equations (77,79) reduce to 

Kl(RtfV) 



M A~^'R A a S ( t ) - _V R^" B 

/" 3 Z A n(r,t) Z A Z B 

I I T? class -J / ; / ; I r> class T> class I 

J l^Aa r l 3=1)3=1 ' Aa ^Bf} I 

and the electronic KS equations simplify to 



(82) 



i 

K N 



9 , x / 1„ 2 i \ f ,3 / n ( r '' t ) 



-EE | r H^»mi +^[".{RgTW}](r,t))yj(r,t) • (83) 
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Equations (82) and (83) are coupled and have to be solved simultaneously. This 
scheme has been applied rather successfully to describe the melting of bulk sodium 
[38]. Compared to the Car-Parrinello method [44, 45, 46] the scheme has the advan- 
tage of not requiring the imposition of orthonormality constraints in the electronic 
equations of motion. 

One has to emphasize that Eqs. (82) and (83) do not involve the Born-Oppen- 
heimer approximation although the nuclear motion is treated classically. This is an 
important advantage over the quantum molecular dynamics approach [47, 48, 49, 50, 
51, 52, 53, 54] where the nuclear Newton equations (82) are solved simultaneously 
with a set of ground-state KS equations at the instantaneous nuclear positions. In 
spite of the obvious numerical advantages one has to keep in mind that the classical 
treatment of nuclear motion is justified only if the probability densities n J 4 Q (R, t) 
remain narrow distributions during the whole process considered. The splitting of 
the nuclear wave packet found, e. g., in pump-probe experiments [55, 56, 57, 58] 
cannot be properly accounted for by treating the nuclear motion classically. In this 
case, one has to face the complete system (67 - 72) of coupled TDKS equations for 
electrons and nuclei. 

4 Electrons in time-dependent electromagnetic fields 

Up to this point we have exclusively dealt with time-dependent electric fields. The 
objective of the present chapter is to incorporate magnetic effects. For simplicity, 
only the electronic degrees of freedom are being discussed, i.e., the nuclear motion 
is not considered. Magnetic fields couple both to the spin and to the electronic 
orbital currents. Hence, the most general TDDFT should encompass both of these 
couplings at the same time. However, to keep matters as simple as possible, we 
shall treat the two couplings separately in the following sections. 

4.1 Coupling to spin 

In order to account for the coupling of a magnetic field B(r, t) to the electronic 
spin, the external potential 

N 

V{t) = Y, «fo . *) = / d3r ^ r Mr, *) (84) 
represented in terms of the density operator 

JV 

n(r)=2«y(r- rj -) (85) 

has to be complemented by a Zeeman term, i.e., V(t) has to be replaced by 

V B {t) = Jd 3 rh(r)v(r,t) -Jd 3 r m(r)-B(r,t) (86) 

where m(r) represents the operator of the spin magnetization. For simplicity we 
assume that the vector B has only one non- vanishing component, the z-component, 
so that 

V B {t)= Jd 3 rn(r)v(r,t)- jd 3 rm z {r)B z (r,t). (87) 
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If the system contains iVj spin-up electrons and N± = N — spin-down electrons 
we can define spin-up and spin-down density operators by 

n T (r) := ^(r-Tj) (88) 
i=i 

JV 

n;(r) := ]T «J(r - r,) . (89) 

j=jV T +l 

In terms of these operators the total density n(r) and the magnetization m z (r) can 
be expressed as 

n(r) = n T (r)+n x (r) (90) 
m z (r) = -/i [n T (r) -n x (r)] (91) 

where /io is the Bohr magneton. 

Defining furthermore the spin-up and spin-down potentials 

Vi(r,t)=v(r,t)+voB z (r,t) (92) 

and 

V l (T,t)=v(T,t)-IM)B x (T,t), (93) 

Eq. (87) simplifies to 



v B (t) 



/A-« t (.M(r,«) + /A-» 1 (.M(r.O. («) 



Starting from the time-dependent many-body Schrodinger equation 

d_ 
dt 



l — y(t)=(f + U + V B (tj)v(t), (95) 



a time-dependent HKS formalism can be established [59] in analogy to section 2: 
The time-dependent spin densities 

n T (r,t) = <*(i)|n T (r)|tt(t)) (96) 
n i (r,t) = <*(t)|n 1 (r)|*(t)) (97) 

evolving from a fixed initial many-body state ^(t ) are in 1-1 correspondence with 
the potentials («f(r, t),u^(r,t)) provided that the latter can be expanded in Taylor 
series around the initial time t . The spin densities thus determine the potentials 
Uf = U|[n|,nJ, «| = «|[n|,nj uniquely up to within purely time-dependent (r- 
independent) additive functions. Consequently the many-body wave function can 
be considered as a functional ^(t) = \P[nf, rojj(i) of the spin densities which is 
unique up to within a purely time-dependent phase factor. Furthermore, following 
the arguments in section 2.2, the spin densities of a given interacting system can be 
determined variationally by solving the Euler-Lagrange equations 

5A B [n hni ] = Q (9g) 



5n T (r,i) 
5A B [n-[,n{\ 



5ni(r,t) 

where the action functional is formally defined as 

rti 



0, (99) 



[n t ,nj := jf 'dt <*[n T , nj(t)|i^ —f—U— Vb|*[n T ,nJ(t)) . (100) 
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Finally, assuming non-interacting v-representability, the spin densities of the inter- 
acting system of interest can be obtained from time-dependent spin orbitals 

N a 

n CT (r,t)=]T|^ CT (r,t)| 2 (101) 
i=i 

coming from the time-dependent KS equations 

*^J>( r '*) = (-^V 2 +v B< > T ,ni](r ) t)^ <Pja(r,t) (102) 

j = l...N„, a =T| 

where the spin-dependent effective single-particle potential for electrons with spin 
a = j, | is given by 

/n(r' t) 
dV— — '— + u xccr [n T ,nJ(r,f). (103) 

Eqs. (101) - (103) constitute the KS scheme of time-dependent spin-density func- 
tional theory. With the xc action functional *4 XC [n T ,nJ defined in analogy to 
Eq. (43), the spin-dependent xc potentials can be represented as functional deriva- 
tives: 

" ro[nT ' nJ(M) = W (104) 

In the limit of vanishing magnetic field the external potentials in Eq. (103) become 
identical 

v^{r,t) =v i {r,t) =v{r,t) for B(r,t) = 0. (105) 

Nevertheless, Eqs. (102) and (103) do not necessarily reduce to the ordinary TDKS 
equations (40) and (41) in this limit. This is because the spin-dependent xc poten- 
tials i> xc j and Vxci are not identical except for the case of spin-saturated systems 
(n T = nj. 

4.2 Coupling to orbital currents 

In order to describe the coupling of time-dependent magnetic fields to the electronic 
orbital currents, the kinetic energy T has to be replaced by 

N | 

(106) 



1 / i \ ' 

^A(t) = ^2 (-<Vr, + -A( rj -,t)J 
.5=1 



where A(r, t) is the time-dependent vector potential related to the magnetic field 

by 

B(r,t) = VxA(r,(). (107) 

Since the vector potential is not a gauge-invariant quantity, particular attention 
has to be paid to gauge transformations: If . . . , rjv,i) is a solution of the 

time-dependent Schrodinger equation 



d 

i 



dt 

then the transformed wave function 

N 



*(*) = (f A (t) + t) - + V'(t)) , (108) 



*(!•!,..., rjv ,t) =exp<j --J2A( Tj ,t) \ V(r u . . . ,r N ,t) (109) 
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is a solution of (108) with the gauge-transformed potentials 3 



v(r,t) = w(r,i) + ~A(r,t) (110) 
A(r,t) = A(r,t) + VA(r,t). (Ill) 
The physical (i.e. gauge-invariant) current density is given by 

j(r,i) - (*(t)|i P (r)|*(t)) + -n{v,t)A{v,t) (112) 

where j p (r) is the paramagnetic current density operator defined in Eq. (14). With 
these preliminaries, the central Hohcnbcrg-Kohn-like theorem [24, 60] to be proven 
subsequently can be formulated as follows: 

The current densities j(r, t) and j'(r, t) evolving from a common initial state 
*o = *(ri, . . . , rjv,io) under the influence of two four-potentials (v(r, i), A(r, t)) 
and (v'(r, t), A'(r, t)) which differ by more than a gauge transformation with A(r, to) = 
are always different provided that the potentials can be expanded in Taylor series 
around the initial time to. 

Since the current density is gauge invariant the proof of the theorem can be 
carried out with an arbitrary representative of the gauge class of (v, A) and an 
arbitrary representative of the gauge class of (v',A'). As representatives we choose 
those four-potentials having a vanishing electric potential, i.e., for v(r,t) we make 
a gauge transformation (110) satisfying 

^A(r,t) = -cv(r,i), A(r,t )=0 (113) 

and for v'(r,t) we make a gauge transformation satisfying 

^A'(r,t) = -ct/(r,t), A'(r,*o) = 0. (114) 

The corresponding gauge-transformed vector potentials are denoted by A(r, t) and 
A'(r,i). Thus we have to show that 

(o,A(r,t)) ± (o,A'(r,t)) (115) 

implies 

j(r,t)^j'(r,f). (116) 

If A(r, to) ^ A'(r,to), then the statement of the theorem is trivially true because 
the initial paramagnetic currents and the initial densities are identical so that 

j(r, t ) - j'(r, t ) = - c n{r, t ) (A(r, t ) - A'(r, t )) + . (117) 

If A(r, t ) = A'(r,t ) the potentials must differ in some higher Taylor coefficient, 



i.e., 



^(A (r ,t)-A' (r ,t)) 1=1 ; *<j (118) 



3 In the context of electrodynamics, the gauge transformation (111) is usually complemented 
by the transformation (p(r,t) = (p(r,t) — l/cdA(r,t)/dt where <p(r, t) is the electric potential. In 
quantum mechanics, on the other hand, one works with the potential energy v(r, i) = q(p(r, t). For 
electrons (q = — e) the gauge transformation of <f> then leads to (110) in atomic units. 
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must be satisfied with an integer I > 0. Calculating the l-th time derivative of the 
densities j(r, t) and j'(r, t) by applying the Heisenberg equation of motion I times 
and taking the difference at the initial time to we obtain 

i|)'(K'.*)-J'(',*))U = ^(',<»)) (M^-A'^^O, 

(119) 

where n(r,t n ) is the particle density at to- By virtue of (119) the current densities 
j(r, t) and j'(r,t) will become different at times infinitesimally later than t . This 
completes the proof. As a consequence of this theorem the physical current density 
j(r, i) determines the potentials t>[j], A[j] uniquely up to within a gauge transfor- 
mation (110), (111). Hence, by virtue of the Schrodinger equation (108), the many- 
body wave function is a current-density functional ^[^(ri, . . . , rjv,i), unique up to 
within a gauge transformation (109). In a fixed gauge, of course, v, A and <]/ are 
determined uniquely by the current density. Applying the theorem to noninteract- 
ing particles then, once again, the potentials v s [j], A s \ji\ and the Slater determinant 
<&[j] leading to the current density j(r, t) are uniquely determined in a fixed gauge. 

In order to derive a TDKS scheme we consider a particular interacting system 
with current density jo (r, t), produced by the external potentials Wo(r, t), Ao(r, t) 
(in a given gauge). Assuming noninteracting w-representability, i.e., assuming the 
existence of potentials v Sy o, A Sj0 leading to jo, we can calculate jo from the equations 

i^j{v,t) = ^-iV+^A, i0 [jo](r,t)^ +«., [jo](r,t)^ J -(r,t) (120) 
1 N 

jo(r,t) = -^(^(r,i)V^ fe (r,t)-^(r,t)V^(r,t)) 

1 k=l 

+ \ (j2\Mr,t)\ 2 ^j A B , (r,t). (121) 

Once the existence of v Si0 and A Si o is assumed, uniqueness follows from the above 
theorem. Up to this point the time-dependent HKS formalism is quite similar to the 
case without magnetic fields developed in section 2. The variational representation 
of v Si o and A Si o, however, turns out to be much more complicated. Following 
Wacker, Kiimmel and Gross [60] the quantum mechanical action functional 

a„,a„[j] - J' 1 dt (mm^ t - *Ao(t) -u- Mwm)) (^) 

has a stationary point at the true current density j , i.e., the latter can be deter- 
mined from the variational equation 



£ A , a Li] 



Sj(r,t) 

Correspondingly, the action functional 



= 0. (123) 

jo 



-4L,A B0 [j] = J* 1 dt (mm§- t - T Asa (t) - v s0 (t)m}(t)) (m 

of noninteracting particles moving in the external potentials v s0 , A s0 has a station- 
ary point at jo as well, i.e., 



Sj(r,t) 



= 0. (125) 

jo 
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In order to deduce an integral equation determining the potentials i; s0 , A s0 , we 
decompose -4u ,A [j] m to a universal part, £>[j], and a functional Q Uo ,A |j] that 
depends on the external potentials vo(r,t), A (r,i): 

X 0l A [j]=B[j]-e wo> Ao[j]- (126) 
The universal part is given by 

= J* 1 dt — T— UMW)) , (127) 

where T is the kinetic energy (6). In terms of the functionals 

n[j](r,t) = - / di'divj(r,i') (128) 

and 

j p [j](r,t) = j(r,i) - l -n[i]{r,t) • A[j](r,i) (129) 
the non-universal contribution to _4[j] can b e expressed as 

Q« ,A„Li] =^ 1 dt J d 3 r |^o(r,i) + ^A (r,t) 2 ^n[j](r,t) + iAo(r,t).j p [j](r,t)| . 

(130) 

Similarly, the action functional of noninteracting particles can be written as 

^ s0 ,A s0 [j]-^[j]-QL,A a0 [j]- (131) 

where 

£ s [j] = J\t{mm^-T\m{t)) 

Ql s0 ,A s0 U] - j t 1 dtj A {(" S o(M) + ^A s0 (r,fj n [j](r,i) 

+ ^A s o(r,t)-j pB [j](r,f)| . (133) 
Note that the functional 

j P .[j] = j(r,t) - -n[i](r,t)A B [j](r,t) (134) 

c 

is different, in general, from the functional j p [j] given by Eq. (129). Defining the 
universal xc functional as 

Ac[j]=£ s [j]-fi[j] (135) 
the total action functional of the interacting system can be expressed as 

A„A [j] - B s []] - Q Vo , Ao [j] - Ac[j] • (136) 

Equating the functional derivatives in (123) and (125) and inserting the expressions 
(131) and (136) we obtain 



$j(r,t) 



jo «('.*) 



(137) 



This equation defines the TDKS potentials i> s o, A s o implicitly in terms of the func- 
tionals A[j] and A s [j]. Clearly, Eq. (137) is rather complicated. The external- 
potential terms Q and Q s are simple functionals of the density and the paramagnetic 
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current density. The complexity of Eq. (137) arises from the fact that the density, 
Eq. (128), and the paramagnetic currents, Eqs. (129), (134), are complicated func- 
tionals of j. Hence a formulation directly in terms of the density and the param- 
agnetic current density would be desirable. For electrons in static electromagnetic 
fields, Vignale and Rasolt [61, 62, 63] have formulated a current-density functional 
theory in terms of the density and the paramagnetic current density which has been 
successfully applied to a variety of systems [63] . A time-dependent HKS formalism 
in terms of the density and the paramagnetic current density, however, has not been 
achieved so far. 

Several extensions of the formalism presented here have been proposed to deal 
with more general situations. Those include superconductors in time-dependent 
electromagnetic fields [60, 64] and time-dependent ensembles either for the electrons 
alone [65, 66] or for the coupled system of electrons and ions [42, 67]. As long as the 
number of photons is large, i.e., > 1 in a volume given by the wave length cubed, 
the electromagnetic fields can be treated as classical fields. For smaller photon 
densities the quantum nature of electromagnetic radiation becomes important. In 
this case, a time-dependent functional theory can be formulated [68] on the basis of 
quantum electrodynamics. In this formulation the electromagnetic field is treated 
as a quantum field to be determined sclf-consistently with the four-current vector 
of the Dirac matter field. 

5 Perturbative regime, basic equations 
5.1 Time-dependent linear density response 

In this section we shall derive [69] a formally exact representation of the linear 
density response ni(rw) of an interacting many-electron system in terms of the 
response function of the corresponding (non-interacting) Kohn-Sham system and a 
frequency-dependent xc kernel. 

We consider electronic systems subject to external potentials of the form 

, ,s / v (r) ; t < t , 1QCA 

V ^ t) = {v (r) +Vl (r,t) ; t > t ^ 

where «o(r) denotes the static external potential of the unperturbed system (typi- 
cally the nuclear Coulomb attraction) and vi(r, t) is a time-dependent perturbation. 
We assume, that at times t < to the system is in the ground state corresponding to 
i>o(r). In this case, the initial density n (r) can be obtained from the self-consistent 
solution of the ordinary ground-state Kohn-Sham equations: 

(-^V 2 + «o(r) + / dV ^7j + ^ n °](r)) &(r) = ^(r) , (139) 

JV 

n (r)=]>>,(r)| 2 . (140) 

By virtue of the static HK-theorem, the initial many-body ground state is uniquely 
determined by the initial ground-state density n . Hence, in this case, the time- 
dependent density n(r,t) is a functional of the external potential alone, 

n(r,t) = n[v ext ](r,t) , (141) 

i.e., there is no additional dependence on the initial many-body state. By virtue 
of the fundamental 1-1 correspondence between time-dependent densities and time- 
dependent potentials, proven by Runge and Gross [7], the functional n[w ex t] can be 
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inverted, i.e., 

Wext(r,t)=«ext[n](r,t). (142) 

Within the realm of perturbation theory, i.e., for sufficiently small vi(r,t) the func- 
tional n[t>ext] can be expanded into a functional Taylor series with respect to the 
perturbation Vi(r, t), 

n(r,t)-n {r)=n 1 {r,t) + n 2 (r,t)+n 3 {r,t) + ... , (143) 

where the lower indices denote the orders in v\ . The first order density response n-i 
is given by 

m(r,i) = Jdt' |dV X (r,t,r',t')Di(r',t') (144) 
with the density-density response function 



x(r>t y )t0 _*»MM) 



5v ext (r',t') 



(145) 



Owing to the static HK theorem, the initial potential v — v cxt [no] is a functional 
of the unperturbed ground-state density n , so that the response function x, by 
Eq. (145), is a functional of n as well. 

For non-interacting particles moving in external potentials v s (r,t), the Runge- 
Gross theorem holds as well. Therefore the functional 

n(r,t)=n[v s ](r,t) (146) 

can be inverted, 

v s (r,t)=v s [n](r,t), (147) 

and the Kohn-Sham response function, i.e., the density-density response function 
of non-interacting particles with unperturbed density no is given by 



/ , / .m Sn[v s ](r,t) 
Xs(r,t,r ,t ) = 



5v s (r>,t>) 



(148) 

f a [no] 



By inserting the functional (141) into the right-hand side of Eq. (147) one has 
formally constructed a unique functional fsbcxt] such that the time-dependent den- 
sity of noninteracting particles moving in v s (r,t) is identical with the density of 
Coulomb- interacting particles moving in u cxt (r,t). The potential v s (r, t) corre- 
sponding to a given w oxt (r,t), is the time-dependent Kohn-Sham potential (41): 

v s (r,t)=v ext (r,t) + Jd 3 r'j^lA +Vxc ( ryt ). (149) 

By virtue of the functional chain rule, the functional derivative of v s with respect to 
u ex t provides a link of the interacting response function (145) to its noninteracting 
counterpart: 



X {r,t,r',t')= Jd 3 x Jdr 



Sn(r 7 t) Sv s (x,t) 



Sv s (x,t) Sv ext (r',t') 



(150) 



Making use of the functional chain rule once more to calculate the functional deriva- 
tive of v s with respect to w oxt one gets 



Sv s (r,t) 



5v ext (r',t') 



= 5{r - r') 5{t - t') 



n 



5{t — t) Sv xc (r,t)\ 5n(x,t) 
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By inserting (151) into (150) and using the definitions (145) and (148) we end up 
with a Dyson-type equation relating the interacting and the noninteracting response 
functions to each other: 

X {v,t,r',t') = Xs (r,t,r',t') + j d 3 x J dr j d 3 x' j dr' Xs (r,t,x,r) 

x ( ^"J^ +/xcN(x,r,x',r')) x(x',r',r',t'), (152) 

where the so-called time-dependent xc kernel 

8v xc [n](r,t) 



f xc [n ](r,t,r',t') := 



5n(r',t>) 



(153) 



is a functional of the initial ground-state density hq. Equations (152) - (153) are 
the central result of our analysis. In previous work, see e.g. Ref. [70], it has been 
common practice to define ,f xc by Eq. (152). The present derivation of Eq. (152) 
from TDDFT shows that / xc , apart from its relation to the response functions x 
and Xs, can also be represented as the functional derivative of the TD xc potential. 
Multiplying Eq. (152) by the perturbing potential vi(r',t') and integrating over 
r' and t' leads to the time-dependent Kohn-Sham equations for the linear density 
response: 

«i(r,t) = Jdt' Jd 3 r' Xs (r,t,r',t')v Stl (r',t'), (154) 
where the effective potential 

v sA (r,t) = «i(r,i)+ J dV + j d 3 r' j dt' / xc [n ](r, t, r', t') m(r', t') (155) 

consists of the external perturbation v\ and the Hartree- and exchange-correlation 
contributions to first order in the perturbing potential v\. We emphasize that 
Eqs. (154) and (155), postulated in previous work [10, 71, 72], constitute an exact 
representation of the linear density response. In other words, the exact linear den- 
sity response n\ (r, t) of an interacting system can be written as the linear density 
response of a noninteracting system to the effective perturbation w Sj i(r,t). Com- 
bining Eqs. (154) and (155) and taking the Fourier transform with respect to time, 
the exact frequency-dependent linear density response is seen to be 

ni(r,w) = J d 3 y X s(r,y;u;)v 1 (y,uj) (156) 
+ j d 3 y J d 3 y' Xs (r,y;uj) ( | y ^ y/ | + /xc["o](y, ")) My' ■ 

The Kohn-Sham response function Xs is readily expressed in terms of the static 
unperturbed Kohn-Sham orbitals (f>k- 

^ ((>j (r)</>*(r)0*(r')0fc(r') 

Here, (fk,fj) are the occupation numbers (0 or 1) of the KS orbitals. The sum- 
mation in (157) ranges over both occupied and unoccupied orbitals, including the 
continuum states. 

In this section we only dealt with the linear response to time-dependent elec- 
tric fields of systems at zero temperature. The corresponding formalism for sys- 
tems at finite temperature in thermal equilibrium was developed in [73, 74]. The 
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currcnt-density-functional response theory for arbitrary time-dependent electromag- 
netic fields has been worked out by Ng [75]. The exchange-correlation kernel / xc , 
given by Eq. (153), comprises all dynamic exchange and correlation effects to linear 
order in the perturbing potential. Depending on physical context, / xc has different 
names: In the theory of the homogeneous electron gas [70, 76, 77, 78] the Fourier 
transform, f xc (q, u)), of / xc (r, t, r', t') with respect to (r — r') and (t — t') is propor- 
tional to the so-called local field correction 

G(q,u) = -^f xc (q,uj). (158) 

In the theory of classical liquids [79], / xc plus the particle-particle interaction is 
known as Ornstein-Zernike function. In practice, of course, this quantity is only 
approximately known. Suitable approximations of / xc will be discussed in section 
6. In order to construct such approximate functionals, it is useful to express f xc 
in terms of the full response function \- An exact representation of f xc is readily 
obtained by solving Eq. (144) for v\ and inserting the result in Eq. (155). Eq. (154) 
then yields 

f xc [n ](r,t,r',t') - x; 1 W(M,r',t')-x- 1 H(r,i,r',t') - , (159) 

where x^ 1 an d X _1 stand for the kernels of the corresponding inverse integral opera- 
tors whose existence on the set of densities specified by Eqs. (138) and (144) follows 
from Eq. (25), as mentioned in section 2.1. The frequency-dependent response oper- 
ators x(r, r'; u) and Xs( r > r '; w), on the other hand, can be non-invertible at isolated 
frequencies [80, 81]. Ng and Singwi [73, 82] have argued, however, that these ex- 
amples are typical of finite systems while for large systems in the thermodynamic 
limit invertibility of the frequency-dependent response operators is guaranteed by 
the second law of thermodynamics. 

As a consequence of causality, the response functions x(r, t, r', t') and x s (r, t, r', t') 
vanish for t' > t. The same statement holds true for the kernels x _1 (r, i, r', t 1 ) and 
XsT 1 ( r J t, r', t 1 ) of the inverse response operators and hence, by Eq. (159), the xc 
kernel must satisfy 

/ xc (r,i,r',i') =0 for t' > t . (160) 

In particular, / xc (r, t, r', t') is not symmetric under exchange of (r,t) and (r',t'). 
Hence, by virtue of Schwarz' lemma for functionals [83], / xc (r, i, r', i') cannot be a 
second functional derivative 5 2 A^ c /Sn(r,t)Sn(r',t'). Since, on the other hand, / xc 
is the functional derivative of v xc one concludes that the exact v xc [n] (r, t) cannot be 
a functional derivative. This is in contradiction to the stationary action principle 
described in section 2.2. which leads to the representation (42) of functional 
derivative. This contradiction is currently an unresolved problem. It appears that 
causality somehow has to be taken into account explicitly in the variational principle 
(36). We emphasize once more that these considerations do not affect the validity 
of the TDKS equations (39) - (41) nor do they affect the validity of the response 
equations (152) - (156). Only the variational representation (42) of w xc appears 
doubtful. 

We finally mention that the chain of arguments leading to Eq. (152) can be 
repeated within static HKS theory as well. This yields 

X stat (r,r') =xf at (r,r') 

+ fd 3 xj dV X f at (r,x) ^^l-^ + /^[n ](x,x'))x stat (x',r'), (161) 



23 



where x stat and xf at are the full and the KS response functions to static perturba- 
tions and 



/ xc [n ](x,x') := 



5v xc [n] (x) 



S 2 E xc [n] 



(162) 



Sn{x') 

On the other hand, taking the Fourier transform of Eq. (152) with respect to (t — t 1 ) 
one obtains 



+ 



d 3 3 



dVxs(r,x;w) 



x — X' 



X(r,r';w) = x s (r,r';w) 
/ xc [n ](x,x';w) ) X(x',r';w) . 



(163) 



Subtracting the zero-frequency limit of this equation from Eq. (161) and using the 
fact that 



X s 



l (r,r') 

W) 



Xs(r, r';w = 0) 
X(r,r';w = 0) 



(164) 
(165) 



one concludes that 



1 d 3 a; y rfV Xs (r,x;a; = 0) (/ xc [n ](x, x'; w - 0) - /^KKx, x')) x(x', r'; w = 0) = . 

(166) 

Acting on this equation with x^ 1 from the left and with x _1 from the right, one 
obtains the rigorous identity [84] 



lim / xc [ra ](x,x';w) 
u)— »o 



S 2 E xc [r 



(167) 



5.2 Time-dependent higher-order response 

Recently there has been a great deal of interest in nonlinear phenomena, both from 
a fundamental point of view, and for the development of new nonlinear optical 
and optoelectronic devices. Even in the optical case, the nonlinearity is usually 
engendered by a solid or molecular medium whose properties are typically deter- 
mined by nonlinear response of an interacting many-electron system. To be able to 
predict these response properties we need an efficient description of exchange and 
correlation phenomena in many-electron systems which are not necessarily near to 
equilibrium. The objective of this chapter is to develop the basic formalism of 
time-dependent nonlinear response within density functional theory, i.e., the calcu- 
lation of the higher-order terms of the functional Taylor expansion Eq. (143). In 
the following this will be done explicitly for the second- and third-order terms 

n 2 {x) = ^fdyfdy'x {2) {x,y,y , )v 1 {y)v 1 {y') (168) 

n 3 (x) = ^[dy Jdy' Jdy" X {3) (x,y,y',y")v 1 (y)v 1 (y')v 1 (y"). (169) 

The extension to terms of arbitrary order is straightforward. For convenience, we 
use the four- vector notation 



x = (r, t) and j dx = j d z r J dt . 



(170) 



The second-and third-order response functions of the interacting system are formally 
given by the functional derivatives 



X {2) (x,y,y') = 



S 2 n(x) 



Sv CKt {y)Sv cxt (y') 
S 3 n(x) 



"0 



fow (y)8v ext (y')Sv ext {y") 



(171) 
(172) 
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of the time-dependent density with respect to the time-dependent external potential 
v oxt evaluated at the initial ground-state density no- From ordinary time-dependent 
perturbation theory, these quantities are given by [85] 

x {2) (x,y,y') = H) 2 £0(*-r)0(i-r') 

V 

(*o|[[nH(x),nH(l/)],nH(l/')]l*o> (173) 

X^{x,y,y',y") = {-if ^ 6{t - r)6(t - T r )6(t' - t") 

v 

(* o I [[[n H (x) , n H (y)} , n H (j/)L % (»")] |*o) (174) 

where the sum has to be taken over all permutations of y,y',y" and the index 
H denotes the Heisenberg picture corresponding to the unperturbed Hamiltonian. 
From the timc-translational invariance of the unperturbed system it follows that 
the response functions (145), (171) and (172) only depend on the differences of the 
time-arguments. Obviously, the full response functions (171) and (172) are very 
hard to calculate. 

The response functions of systems of noninteracting particles, on the other hand, 
are functional derivatives of the density with respect to the time-dependent single- 
particle potential v s : 



(2) , S 2 n(x) 

5v s (y)Sv s (y') 

y {3) (x v v v") = <53n(x) 

Xs [x,y,y,y ) 6 v s (y)Sv s (y')Sv s (y") 



(175) 
(176) 



These functions can be expressed in terms of single-particle orbitals, similar to the 
linear response function (157). 

To obtain the higher-order expressions of the density response, we use the func- 
tional chain rule in Eq. (171): 

5 2 n{x) 5 f 5n(x) Sv s (z) 



Sv cxt (y)Sv cxt (y') 5v cxt (y) J 5v s (z) Sv cxt (y') 

5 2 n(x) Sv s (z') Sv s (z) 



= j dz J dz' 



8v s (z')5v s (z) Sv ext (y) Sv ext (y') 



dJj^l gVg) (177) 

Sv s (z) Sv cxt (y)v cxt {y') 

As has been outlined in Section 2, the full time-dependent Kohn-Sham potential 

v B (x) = v cxt (x) + v H (x) + v xc (x) (178) 

is a unique functional of the external potential v cxt . Hence, we get 

S 2 v s (z) S /" , S(v H (z) +v xc (z)) 5n(z') 



I 



Sv cxt (y)Sv cxt (y') 6v cxt (y) J " Sn(z') Sv cxt (y') 

J 2 (v H (z) + v xc (z)) Sn(z") Sn(z') 



dz' / dz 



5n(z")Sn(z') 6v cxt (y) Sv cxt (y') 

J Sn(z') Sv ext (y)5v ext (y') ' 

Combining Eqs. (177) and (179) and evaluating all functionals at the ground-state 
density no, we obtain a Dyson- type relation for the second-order response function 
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(171): 



x (2) {x, y, y 1 ) = JdzJ dz' x f ] (x, z,z') 



Svs(z') 
no Sv eKt (y') 



+ J dz Xs (x, z) J dz' J dz" g xc (z, z', z") X (z', y) X {z" , y') 

+ I dz Xs (x,z) ( dz' (w(z,z') + U c (z,z'))x i2) (z',y,y') ,(180) 



(181) 



where the time-dependent second-order xc kernel g xc is defined as: 

/ / //n S 2 v xc (z) 
g xc (z,z ,z ) = 



5n(z')Sn(z") 



To arrive at Eq. (180) we have used the definitions (145), (148), (171) and (175) 
of the density response functions. Furthermore, we have abbreviated the kernel of 
the (instantaneous) Coulomb interaction by w(x,x') = 6(t — t')/\r — r'\. Finally, 
by inserting Eq. (180) into (168) one arrives at the time-dependent Kohn-Sham 
equations for the second-order density response: 



n 2 ( 



(x) = - dz dz' xi 2) (^^> s ,i(z)v M (z') 



+ 



\ j dz j dz' j dz" Xs(x,z)g xc (z,z',z")ni(z')ni(z") 



+ dz dz' x s (x,z)(w(z,z') + Uc{z,z'))n 2 {z') . 



(182) 

Solving Eqs. (154) and (155) first, allows for the subsequent solution of the selfcon- 
sistcnt Eq. (182) which is quadratic in the (effective) perturbing potential (155). 

In similar fashion, one can set up the equation for the third-order density re- 
sponse (169): 

Mx) - \ j dy j dy' J dy" X f ] (x, y, y' , y")v sA (y)v sA (y')v sA (y") 

2~ J dV j dy ' j ' dZ J dZ ' x « 2) ^' y ' y '^ s ' 1 ^ 5xc ^'' z ' z ') ni ( z ) ni ( z ') 
dy J dy' J dy" X ?\x,y,y')v s ,i{y) (w(y',y") + f xc (y',y"))n 2 (y") 



+ 
+ 



dy dz dz 



z' J dz" Xs{x, y)h xc (y, z, z' , z")n 1 (z)n 1 (z')n 1 {z") 



dy J dy' J dy" x s {x,y)g xc (y,y' ,y")ni(y')n 2 (y") 
dy / dy' x s {x,y) {w(y,y') + f xc (y,y'))n 3 (y') . 



(183) 



The quantity h xc occurring in this equation is the third-order functional derivative 
of the time-dependent xc potential with respect to the time-dependent densities: 



K c {y,z,z' ,z") 



S 3 v xc (y) 



Sn(z)Sn(z')Sn(z") 



(184) 



Interestingly, the equations (154), (155), (182) and (183) for the z-th order density 
responses all exhibit the same structure: 

rii{x) = M l {x) + J dy J dy' x s {x,y)(w(y,y') + f xc (y,y'))rii(y') i = 1, 2,3 , 

(185) 
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where the functionals A4i(x) are known after the solution of the (i — l)th order. 
This establishes a hierarchy of Kohn-Sham equations for the time-dependent density 
response. 

The frequency-dependent nonlinear density responses are given by the Fourier 
transforms of Eqs.(185). For monochromatic perturbations, the expressions for the 
higher-order frequency dependent density shifts decouple in the frequency variable. 
The corresponding formulae and explicit expressions for the Kohn-Sham response 
functions up to third order are given in work of Senatore and Subbaswamy [86]. 
The corresponding static higher-order response has been worked out and applied to 
solids by Gonze and Vigncron [87]. 

6 The time-dependent exchange-correlation poten- 
tial: Rigorous properties and approximate func- 
tionals 

6.1 Approximations based on the homogeneous electron gas 

The simplest possible approximation of the time-dependent xc potential is the so- 
called time- dependent or "adiabatic" local density approximation (ALDA). It em- 
ploys the functional form of the static LDA with a time-dependent density: 

v^ A [n](r,t) = « x h c om (n(M)) = A (pe x ° m (p)) | p= „ (r , t) ■ (186) 

Here e x J? m is the xc energy per particle of the homogeneous electron gas. By its 
very definition, the ALDA can be expected to be a good approximation only for 
nearly homogeneous densities, i.e., for functions n(r,t) that are slowly varying both 
spatially and temporally. It will turn out, however, that the ALDA gives rather 
accurate results even for rapidly varying densities (see sections 7 and 8). For the 
time-dependent xc kernel (153), Eq. (186) leads to 

/ x A c LDA [n ](r,t,r',t') = S(t - t')5(r - r')^ (p4r(p))\ p=no(r) • (187) 
The Fourier-transformed quantity 

/ x A c LDA [n ](r,r';c) = 5(r-r')^ (pe^(p)) U o(r) • (188) 

has no frequency-dependence at all. 

In order to incorporate the frequency-dependence of / xc in some approximate 
fashion, Gross and Kohn [71] suggested to use the frequency-dependent xc kernel 
/ x ° m of the homogeneous electron gas in the sense of an LDA: 

/ x L c DA [n ](r,r';c) := / x h c om (n (r), |r - r'|; w ) . (189) 

The LDA of non-local quantities, such as response functions, always involves 
some ambiguity [1, 2] as to whether the inhomogeneous n is to be evaluated at 
r, at r', or at some suitably chosen mean value of r and r'. Of course, in the 
limit of slowly varying no(r) (i. e. in the limit where the LDA should be a good 
approximation) the choice does not matter. In addition to the LDA replacement 
/xc - * /xc " 1 ! Gross and Kohn [71] made the assumption that ni(r, uj) is slowly 
varying on the length scale given by the range of / xc om (no (r) , |r — r'| ; w) . Under this 
assumption, the change in the xc potential can be calculated as 

v£{T,u)= ni {r,oj) jd s r' /£">o(r),|r-r'|;u;) . (190) 
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In terms of the Fourier transform of f^° m with respect to (r— r'), Eq. (190) amounts 
to the approximation 

f^[n }(v,v';u;) = S(r-v')^r(n (v),q = 0;uj) . (191) 

This approximation requires the xc kernel of the homogeneous electron gas as input. 
In order to investigate this quantity we consider Eq. (159) in the homogeneous case. 
Fourier transformation with respect to (r — r') and (t — t') leads to 

/4 om (n , q; w) - — h 7^ r - ~ h ^ r - % • (192) 

The response function Xs° m °f a non-interacting homogeneous system is the well- 
known Lindhard function. The full response function x hom \ on the other hand, is 
not known analytically. However, some exact features of X hom are known. From 
these, the following exact properties of f^° m can be deduced: 

1. As a consequence of the compressibility sum rule one finds [76] 

lim tfr(q, w = 0) = f- 2 {ne h x ° m {n)) = /o(n) . (193) 

This shows that /£ LDA , as given by Eq. (188), is identical with the zero- 
frequency limit of /^ K . 

2. The third-frequency-moment sum rule leads to [88] 

lim.£ m (^ = «)) 

3. According to the best estimates [89, 90] of ej™, the following relation holds 
for all densities: 

/o(n) < /oo(n) < . (195) 



4. The large-g behavior at zero frequency is given by [91] 

lim fZTfau, = 0) = -2n^ (V 1 ^^)) % 
q^oc an \ / <? z 



The function B(n) has been fitted [92] to Monte-Carlo results. The resulting 
parametrization 

1 + 2.15a; + 0.435a; 3 ^_ / 3 \ 1/6 



3- l^-.OTOO,- X = ^ = [^n) 

reproduces the Monte-Carlo results with a precision of about 1% in the density 
range < r s < 10. 

5. The short- wavelength behavior in the high-frequency limit is given by [93, 94] 

lim f%r{q,u> = oo) = -l- %{l , 9 (0)) + 6^ A. ( n -W (i 98 ) 
q^oo 3 an V / 

where g(r) denotes the pair correlation function. 
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6. f^° m (q,uj) is a complex- valued function satisfying the symmetry relations 

$trtr(<i,u)=MfZr(<i,-") (199) 
^r(i^) = -^r(Q,-^) • (200) 

7. f£° m (q,oj) is an analytic function of ui in the upper half of the complex em- 
plane and approaches a real function foo(q) for uj — > oo [70]. Therefore, the 
function (f^° m (q,u}) — foo(q)) satisfies standard Kramers-Kronig relations: 

m h -{^) /cow = p/ v ^r (201) 

^r(^) = -P/y^ 7 "^" • (202) 

8. The imaginary part of /x° m exhibits the high-frequency behavior 

limS/^fe,)^^ (203) 

for any g < 00 [95]. A second-order perturbation expansion [95, 96] of the 
irreducible polarization propagator leads to the high-density limit 

"f ■ « 

Other authors [97, 98] find c = 46tt/15; see also Ref. [94]. 

9. In the same limit, the real part of f^° m behaves like [71] 

lim (q,u>)=f 00 (q) + ^ . (205) 

Since c > 0, the infinite-frequency value is approached from above. This 
implies, in view of the relation (195), that ^H.f^° m (q — 0, uj) cannot grow 
monotonically from / to f^. 

The above features of j^° m are valid for a three-dimensional electron gas. Analogous 
results have been obtained for the two-dimensional case [95, 99, 100]. 

Taking into account the exact high- and low-frequency limits, Gross and Kohn 
[71] proposed the following parametrization for the imaginary part of f^° m : 

* 4 = -^(i^F ■ (206) 

where 

a(n) = - c ( 7 / c ) 5 /3 (/oo(n ) - ioM) 5/3 (207) 

6(n) = ( 7 /c) 4 / 3 (/ oo (n)-/ (n)) 4 / 3 (208) 

..-^m? . (209) 



4V2tt 

fo, /oo, and c are given by Eqs. (193), (194), and (204), respectively. Using the 
Kramers-Kronig relation (201), the real part can be expressed as 

- f 

— JCOT nil , 

7TS Z V 



2£ |i V^n^ 8 1 



l^s n /l + s 1 



2 '^2 



2 J 2 V 2 '72, 

s 2 = 1 + bcu 2 . (210) 
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0) -10 



r s = 2 
r. = 4 



0) (a.U.) 



Figure 1: Real part of the parametrization for / xc om (<7 = 0, w), from Rcf. [88] 



£7 and II are complete elliptic integrals of the second and third kind in the standard 
notation of Byrd and Friedman [101]. This completes the explicit form of the Gross- 
Kohn approximation (191). 

Figs. 1 and 2 show the real and imaginary part of / xc om as calculated from (206) 
and (210). The functions are plotted for the two density values corresponding to 
r s = 2 and r s = 4. For the lower density value (r s = 4), a considerable frequency 
dependence is found. The dependence on to becomes less pronounced for higher 
densities. In the extreme high-density limit, the difference between f and 
tends to zero. One finds the exact result 

foo~fo~r 2 s farr a -0 . (211) 

At the same time, the depth of the minimum of S/^"" 1 decreases, within the 
parametrization (206) proportional to r 2 s . 

We finally mention that an extension of the parametrization (206) to non- 
vanishing q was given by Dabrowski [102]. The spin-dependent case was treated 
by Liu [103]. A similar interpolation for the exchange-correlation kernel of the 
2-dimcnsional electron gas has been derived by Holas and Singwi [95]. 

In the construction and improvement of static ground-state density functionals, 
various exact constraints such as xc hole normalization [104] and scaling relations 
[105] have been extremely useful. While the development of explicit time-dependent 
functionals is at a comparatively early stage, there are some constraint conditions 
which can be useful in the time-dependent context. First of all, some of the ex- 
act properties of the homogeneous-electron-gas kernel /J™ are readily general- 
ized to the inhomogeneous case: Causality leads to Kramers-Kronig relations for 
/ xc (r, r';w) analogous to Eqs. (201) and (202), and the fact that / xc (r, t, r', i') is a 
real- valued quantity implies that 

/ xc (r,r';u;) = / xc (r,r';- W )*. (212) 
Besides that, the response functions \s an d X satisfy the symmetry relations [106] 

X (r,r';uO = x(r',r;w) (213) 
X B (r,r';w) = x B (r',r;w) (214) 

provided that the unperturbed system has time-reversal symmetry. Equation (163) 
then implies that 

/ xc (r,r» = / xc (r',r;o;). (215) 
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o) (a.u.) 



Figure 2: Imaginary part of the paramctrization for f x ° m (q — 0,oj), from Ref. [88]. 

Further exact constraints can be deduced from the quantum mechanical equation 
of motion (15). For the operator 

f = J d 3 rrh(r) (216) 

Eq. (15) leads to 

;<*(t)|r|tt(t)> = jj d 3 rrn(r,t)= l (^(t)\[H(t),m(t)} (217) 



d 

Jt' 



where 



H(t)=f + U + V cxt (t). 



(218) 



Taking the time derivative of Eq. (217) and employing the equation of motion (15) 
once more one obtains 



d^ 
dt 2 



J Arn(r,t) = -<tt(t)| [H(t), [H(t), f]]|*(i)) (219) 



because J^[i?(t),r] = 0. Using the translational invariance of the Coulomb in- 
teraction U the double commutator in (219) is readily calculated, leading to the 
traditional Ehrenfest theorem: 



dt 2 



J d 3 rrn(r,t) = - J d 3 rn(r, t)Vv eKt (r, t) . (220) 

Likewise, for noninteracting systems described by Hamiltonians of the form 

H s {t)=f + V s (t) (221) 

one obtains 



Arn s (r,t) = - / d 3 rn s (r, i)Vv s (r, t) . 



dt 2 

For the unique KS potential 

v s [n] (r , t) = v cxt (r,t) + v H [n] (r, t) + v xc [n] (r , t) 



(222) 



(223) 
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which reproduces the density n(r, t) of the interacting system, Eq. (222) leads to 

^ J d 3 rrn{v,t) = -J d 3 rn(r, t)V (v cxt (r, t) + v H [n](r, i) + u xc [n](r, t)) . 

(224) 

Subtracting Eq. (220) from Eq. (224) one obtains the rigorous result 

J d 3 rn(r, t)Vv xc [n] (r, t)= 0. (225) 
To arrive at Eq. (225) we have used the fact that the Hartree potential 

v H [n](r,t)= /rfVj^ (226) 

satisfies the equation 

J d 3 rn(r, t)Vv H [n] (r, t) = . (227) 

Equation (225) was first obtained by Vignale [107] from invariance properties of the 
xc action functional A xc defined in Eq. (43). The derivation given here [108, 109] 
has the advantage of being independent of the stationary action principle. 

Applying the equation of motion (15) to the angle operator ip and using the 
rotational invariance of the Coulomb interaction U, one obtains 



dt 2 



(*(t)|0|*(t)) =- J d 3 rn(r, t)v x V« cxt (r, t) . (228) 



Subtraction of the corresponding equation for the KS potential (223) then leads to 
the exact constraint 

J d 3 rn(r,t)r x Vv xc [n](r, i) = 0. (229) 

Corresponding properties of the exact xc kernel are obtained by evaluating the 
left-hand sides of Eqs. (225) and (229) at the density 

n(r,t)=no(r) + <fn(r,t), (230) 

where 5n(r,t) is an arbitrary deviation from the ground-state density no(r). To 
first order in 5n one obtains from Eq. (225) 



= J d 3 rn (r)\7vxc[n }(r) 
+ J d 3 r' j dt'Sn{r',t') 



5(t-t')V r >v xc [n }(r') 



+ 



J d 3 rn (r)V r f xc [n }(r,t,r',t') 



(231) 



The first integral on the right-hand side of this equation must vanish. (This is the 
static limit [110] of Eq. (225).) Since Sn is arbitrary, the second integral leads to 
the identity 

J d 3 rn (r)V r / xc [n ](r,t,r',t') = -6{t - t')V r -D xc [n ](r') . (232) 

Taking the Fourier- Transform of this equation with respect to (t — t') one obtains 
the constraint 

J d 3 rn (r)V r / xc [n ](r,r';w) = -V r ^ xc [n ](r') . (233) 
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Applying the same procedure to Eq. (229) one arrives at [108, 109] 

J d 3 rn (r)r x V r / xc [n ](r, r'; u) = -r' x V r ^ xc [n ](r') . (234) 
Finally, multiplying Eqs. (233) and (234) by no(r') and integrating over r' leads to 
J d 3 r J d 3 r'no(r)no(r')V r / xc [n ](r,r'; W ) = (235) 

and 

f d 3 r [ dVn (r)n (r> x V r / xc [n ](r, r'; w) = . (236) 



Equation (233) was first obtained by Vignale [111] from a new sum rule for the 
response function. The ALDA satisfies the constraints (233) and (234) while the 
Gross-Kohn approximation (191) is easily seen to violate them. This fact is closely 
related to the violation of the Harmonic Potential Theorem which will be discussed 
in detail below. 

Another type of constraint on theories of time-dependent phenomena in inter- 
acting inhomogeneous systems is obtained by taking expectation values of repeated 
commutators of current operators with the Hamiltonian. In this way one obtains ex- 
act relations for frequency moments of response functions. Very recently Sturm [94] 
has given a detailed study of the odd frequency moments of the dielectric function 
in inhomogeneous systems, and has explored odd moments up to the seventh. The 
ALDA satisfies the first, third and fifth frequency moment sum rules but violates 
the seventh as was demonstrated by Sturm for metals in the nearly free electron 
approximation. 

Finally, another rigorous constraint [112] is known as the Harmonic Potential 
Theorem (HPT), as it relates to the motion of interacting many-electron systems in 
an externally- imposed harmonic oscillator potential v(r) = |rKr plus the potential 
— F(t)-r describing a spatially uniform, time-dependent external force F(t). Here 
K is a spring-constant matrix which can be assumed symmetric without loss of 
generality. (Suitable choices of K yield various physical situations: for example, 
the choice K = diag(fc, k, k) corresponds to a spherical quantum dot or "Hooke's 
atom", while the choice K = diag(0, 0, k) yields a parabolic quantum well such 
as may be grown in the Gai_ x ALAs system by molecular beam epitaxy.) The 
harmonic external potential is special, being the only confining potential which 
retains its form when one transforms to a homogeneously accelerated reference 
frame. To see that it does so, consider [113] a moving frame whose origin has the 
space coordinate X(i) relative to the rest frame. The observer in this frame sees a 
total external potential 

v(f,t) = ^r- K-r+mX-f - F(i) -r (237) 

where r = r X(i) is the position coordinate in the new frame, and the second 
term in (237) is the centrifugal or fictitious potential due to motion of the frame. 
If X(t) satisfies the classical equation of motion 

mX(t) = -K • X(t) + F(t) (238) 

then one obtains the potential in the moving frame as 

v(r,t) = ^r-K-r + c(t). (239) 

This transformed external potential (239) has the same form as the potential for 
the undriven (F = 0) harmonic well problem in the rest frame, except for the term c 
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which depends on time but not on r. Furthermore, because the (Coulomb or other) 
particle-particle interaction is a function only of differences r, — rj = ?i — fj, the 
interaction potential is also invariant under the transformation to the accelerated 
frame. Thus, both classically and quantum mechanically, any state or motion in the 
rest frame has a counterpart motion with superimposed translation X(t), provided 
that (238) is satisfied. In particular 

for harmonically- confined interacting systems there exist quantum states in which 
the ground-state many-body wavefunction is translated rigidly (up to a phase factor) 
as in classical motion, and hence the ground-state number density no(r) is replaced 
by the rigidly moving density n(r,t) = no(r — X(t)). 

(240) 

This conclusion is the Harmonic Potential Theorem (HPT). It is an extension of 
the generalized Kohn Theorem [114]: the latter only refers to the frequency de- 
pendence of linear response and does not address the spatial profile of the moving 
density. The HPT can be also proved more formally [112] by explicit construction 
of the moving many-body wavefunction as seen in the rest frame. It is important to 
note that systems confined by a scalar harmonic potential (e.g. quantum dots) are 
spatially finite in at least one dimension, and have strong spatial inhomogeneity at 
their edges. Thus the HPT constitutes an exact result beyond the level of linear re- 
sponse for the time- dependent behavior of an inhomogeneous, interacting many-body 
system. As such, it poses an interesting constraint on approximate general theories 
of time-dependent many-body physics, such as local-density versions of TDDFT. Of 
course, since the HPT is valid generally it is also valid for linear response. Vignale 
[107] has shown that the HPT result holds even with the inclusion of a homogeneous 
magnetic field. 

Another closely related constraint is that of Galilcian invariance. Suppose that a 
many-body wavefunction ^(ri, r^, rjv) satisfies the time-independent interacting 
N-particle Schrodinger equation with an external one-particle potential v(r). Then, 
provided that the inter-particle interaction depends on coordinate differences only, 
it is readily verified that a boosted wavefunction of the form 

exp(-iS(t) + iU- rj)*(ri - Ut, r 2 - Ut, rjv - Ut) , (241) 

where S(t) is the corresponding classical action [112], satisfies the time-dependent 
interacting N-body Schrodinger equation with boosted external potential v(r — Ut). 
Because the phase factor disappears in forming , this result implies that all 
many-body probability densities are rigidly boosted when the external one-body 
potential is boosted. In particular, the one-particle density is rigidly boosted, and 
this particular aspect of Galileian invariance should apply to TDDFT which deals 
directly with densities. In applying this criterion, it will clearly be necessary to 
relax the condition used in section 5 that the initial many-body wavefunction be 
the ground-state wavefunction. In fact, to represent a system boosted to constant 
velocity U, the initial wavefunction must contain the additional phase factor shown 
in Eq. (241) 

The following question now arises: which approximations in TDDFT satisfy 
the HPT and Galileian invariance? By noting that the ALDA xc potential rigidly 
follows the density when the latter is rigidly moved, and by examining the TDKS 
equations for harmonic confinement with and without a driving field Dobson [112] 
showed that the ALDA satisfies the HPT for motion of arbitrary amplitude. The 
same proof in fact shows that any approximation to TDDFT satisfies the HPT 
provided that the xc potential rigidly follows a rigidly translated density. This 
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rigid-following condition will be termed Generalized Translational Invariance and 
can be expressed as: 

v xc [n'](r,t)=v xc [n](r-X(t)) (242) 

Here n(r) is an arbitrary time-independent density and n'(r, i) = n(r — X(t)) is 
the same density rigidly boosted. (The otherwise-arbitrary displacement function 
X(t) will need to be zero at the initial time to, and the initial many body state 
will need to be the ground-state, in order for v xc in (242) to be defined in the same 
manner used earlier in this chapter.) Eq. (242) simply says that a rigid (possibly 
accelerated) motion of the density implies a similar rigid motion of the xc potential. 
Equation (242) was first demonstrated by Vignale [107] from the covariance of the 
time-dependent Schrodingcr equation under transformation to an accelerated rest 
frame. Vignale also generelalized the treatment to include a magnetic field. The 
same condition (242) with X(t) = XJt will also ensure that an approximation to v xc 
satisfies Galileian invariance. 

Perhaps surprisingly, the Gross-Kohn approximation (191) unlike the ALDA, 
does not satisfy the HPT constraint. This was proved in Ref. [112] by exhibiting a 
specific counterexample. 

The question now arises how one might correct this situation. One attempt 
[112] is based on the heuristic picture that, in the rigid HPT motion, all the relative 
particle motions [115], and therefore the exchange and correlation phenomena, are 
exactly as in the ground-state. In particular, the static exchange-correlation kernel 
/ xc (n, u = 0) is appropriate for this very special motion, even though the frequency 
of the HPT motion has the high value u = up. This is why the ALDA succeeds 
with the HPT motion: it uses f xc (n 7 uj = 0) in all circumstances and therefore is 
fortuitously exact for HPT motion. The original GK formalism requires the use 
of f xc (n,uj = uip), and this is the core of the difficulty. (A similar difficulty was 
also demonstrated [112] for hydrodynamic theory of plasmons where, once again, 
the frequency dependence of a coefficient, fi 2 (w — ► oo) ^ /? 2 (u> — > 0), is to blame, 3 
being the pressure or diffusion coefficient.) 

While most motions are not simple rigid displacements, there will be an element 
of this type of motion, as well as an element of compression, in more general motions 
provided that the original density is spatially inhomogeneous. In the linear response 
regime, a well-defined separation between these two components of the motion can 
be made [112] by first introducing the fluid velocity u(r, t) — J(r, t) /n(r, t) where the 
exact current density J can be obtained from the TDKS orbitals as demonstrated 
in section 2.3. A fluid element displacement X(r,t)is then defined for a general 
motion by 

u(r, t')dt' , u=— (243) 

and by integrating the linearized continuity equation with respect to t at fixed r 
we obtain an expression for the perturbation to the equilibrium density no(r) in a 
small motion: 

m(r,t) = -V • [no(r)X(r,t)] = -n (r)V • X(r,i) - X(r,t) • Vn (r) . (244) 

For an arbitrary motion, we denote the first term in (244) as 

ni„(r,t) = -no(r)V-X(r,t) (245) 

and interpret m as the compressive part of the density perturbation, to be asso- 
ciated with /xc(w) where w is the actual frequency of the linear motion. The other 
density component from (244) is 

n lb (r, t) = -X(r, t)-Vn (r) (246) 
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and this is the part one would have obtained if the equilibrium density had been 
rigidly translated, suggesting that it should be associated with a zero-frequency 
kernel / xc (w = 0). These two components make up the total density perturbation, 

m = n la + nib , (247) 

and the above arguments suggest that the xc potential for small-amplitude motion 
at frequency to should be 

vi xc {r,u) =/ xc ("o(r),u;)ni (r,a;) + / xc (no(r),w = 0)ni 6 (r,w). (248) 

It is immediately apparent that (248) will give the correct zero-frequency xc poten- 
tial value for Harmonic Potential Theorem motion. For this motion, the gas moves 
rigidly implying X is independent of r so that the compressive part, n la , of the 
density perturbation from (245) is zero. Equally, for perturbations to a uniform 
electron gas, Vn and hence nib is zero, so that (248) gives the uniform-gas xc 
kernel / xc (w) at the actual frequency u), as required. 

A modification similar to (248) was also proposed in [112] for the pressure or 
diffusion term in hydrodynamics, and the resulting formalism has had some success 
with a unified description of boundary conditions and plasmon modes on parabolic 
wells [116]. 

Since the fluid displacement during linear response at a definite frequency ui is 
given by X = J /(icon), the postulated Eq. (248) suggests that v xc is not a local 
function of the density but rather of the current density J. There are, however, pre- 
liminary indications that, for nonlinear phenomena such that a definite frequency 
cannot be assigned to the motion the fluid displacement X may yield a more di- 
rect formulation of xc phenomena than does the current density (see later in this 
chapter) . 

Numerical applications of the new formalism implied by (248) are under develop- 
ment [117]. Preliminary indications are that the ALDA, the Gross-Kohn approxima- 
tion (191) and (248) will all give substantially different results for at least one of the 
plasmon modes of a low-density parabolic quantum well, say for r s = 6. (The modes 
in question are the HPT ( "Kohn" or "sloshing" ) mode, the standing plasmon modes 
[118], and also the 2D plasmon mode at substantial surface-directed wavenumber 
<7|| for which case the frequency is not constrained by model-independent theorems 
[119]-) 

Vignale [107] has given an alternative method to ensure that any xc formalism 
with finite memory (i.e. with frequency-dependent xc kernel) will satisfy the HPT. 
Starting from a simple Ansatz for the action integral, he derived an xc potential 

vixc(r,t)= / f xc (n (r),t-t')Sn re i(r,t')dt' (249) 

where 

5n rel (r,t) = n(r + TL cm (t)) - n (r) (250) 

is the density perturbation seen by an observer moving with the global center of 
mass 

Rcm(t) = ^ J rn(r,t)d 3 r. (251) 

This approach does ensure satisfaction of the HPT. It differs from the method 
described above in that it is very much less local, requiring a determination of the 
global center of mass from (251) at each instant t. One can imagine situations where 
the two formalisms will give very different results. For example, consider two well- 
separated layers of electron gas confined in parallel parabolic wells. At the Hartree 
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and Hartree-Fock levels there is no interaction between these wells in the absence 
of significant wavefunction overlap (and in the absence of any perturbation which 
might break the symmetry in the plane of the electron gas layers). For sufficient 
separations any residual van der Waals interaction can be made as small as desired, 
so the <7|| = modes of oscillation of the two wells will be uncoupled. First consider 
a motion in which the two electron gases execute HPT motion (sloshing sideways) 
in phase. Then the global electronic center of mass also executes HPT motion 
and the Vignale method will give the correct HPT motion of the combined system. 
Secondly, however, consider the mode in which the two sloshing motions are 180 
degrees out of phase. Then the global center of mass is stationary and the Vignale 
correction makes no difference, leaving the GK formalism unmodified. But this 
method is known not give the HPT motion correctly, as discussed above. The 
method described by Eq. (248), on the other hand, is more local in its effect and it 
corrects the motion of each well separately, giving the correct HPT motion of each 
gas even for the 180 degrees phase mode. 

In general, for systems far from equilibrium it is not at all clear how one should 
approximate the full xc potential v xc [n](r 7 t). The most general possible nonlinear 
dependence of v xc [n] (r, t) on n must involve at least terms with n evaluated at one 
space-time point £' = (r' ,t'), terms with n evaluated at two spacetime points £' and 

terms with n evaluated at three points and so on. (Even this might 

not cover all possibilities, but the only counterexamples so far noted have involved 
essentially singular functions [120].) Thus in general we expect to require nonlinear 
functions such that 

«xc[n](o = / dewww),t,?)+ 1 «w( 2 V(eV(r'),«'e")+- (252) 

The functional derivative of (252) is 



(n 



5n(£') dn 

(253) 

where w[ 2) = dW(n',n",Z,£',Z")/dn' and W^ 2) = dWM(n',n",£,£',£")/dn". 
Considerable simplification is achieved by postulating the following local-density 
Ansatz (254) [121] for the functional derivative 

Sv xc (r,t) _ hom 



6n(r',t') 



f*r(n(r',t'),\r-r'\,t-t') (254) 



where /^° m (n, |r — r'| , t — t') is the nonlocal, delayed xc kernel of the uniform elec- 
tron gas of density n. Clearly, in the limit of weakly inhomogeneous systems, i.e., 
for systems with densities n(r',t') — ► const this Ansatz becomes exact. We now 
seek a functional v xc [n](r,t) whose functional derivative Sv xc (r,t)/Sn(r' ,t') is given 
by (254). The task of finding such a v xc is possible because (254) is a function of 
density at one point only. Hence the integral terms in (252) and (253) involving 
and higher must be discarded and it follows that 

^-(n(r', t'),vt, v't') = / x h c om (n(r', t'), |r - r'| , t - t'). (255) 

Thus is a density integral of f x ° m (n, r, t). Assuming that v xc is zero in a 

zero-density system, and defining 

W xc (n, r, r) - f / x h c om (p, r, r)dp (256) 



o 
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wc then obtain 



v KC [n](r,t) = j ' dt'd 3 r'W xc (n(r',t'),\r-r'\,t-t'). (257) 

The Ansatz (257) makes v KC (r,t) depend principally on the density near the 
point r, at a range of times t' which are near to, but earlier than, t. In the following 
we propose to improve this by noting that, if there is streaming in the many-body 
fluid, the memory of past densities is likely to be greatest when one remains with 
the same fluid clement rather than remaining with the same spatial point r. Thus 
we propose [122] instead of (257) 

v xc [n]{r,t) = J dt'd 3 r'W xc (n(r',t'),\K(t' | r,t) -r'| ,t-t'). (258) 

where W xc is still given by (256). 

In (258) the density from a past time t' which most strongly influences v xc (r,t) 
is the density at position r' = R(i' | r,i), where the trajectory function R(t'|r, i) 
of a fluid element is its position at time t' , given that its position at time t is r. 
R can be defined unambiguously by demanding that its time derivative is the fluid 
velocity u = J/n formed from the current density J(r, t): 

^R(i'|r,i) = u(R,t / ) = J(R,t / )/n(R,t / ) (259) 

where all occurrences of R have the same arguments as on the left-hand side of Eq. 
(259). The boundary condition on (259) is 

R(t|r,f) = r. (260) 

In (258) one acknowledges that the physics of delayed correlation will have its 
maximum degree of spatial locality if the observer is riding on a fluid element rather 
than observing from a fixed location r. Eqs. (258) - (260) represent our general 
expression for the dynamic xc potential. The use of (258) in place of (257) will turn 
out to provide a nonlinear theory which, regardless of its validity in other respects, 
at least satisfies both the nonlinear Harmonic Potential Theorem [112] and the 
requirements of Galileian invariancc [107]. To demonstrate that (258) satisfies the 
HPT, we show that it satisfies the generalized Galileian invariance condition (242). 
The only difficulty is that (258) has an implicit and highly nonlocal dependence on 
n(r',t) via the current density dependence of R. From [112], however, it follows 
that for HPT motion the TDKS equations involve not only a rigidly boosted density 
n'(r,t) — n (r — X(t)), but also a boosted current J'(r, t) — n'(r, t)~K(t) because 
of the phase factor introduced by the motion into the KS wavefunctions. Thus the 
fluid velocity is just u(r, t) = X(t). From (259) and (260) it then follows that 

R[n'] {t 1 | r, t) = r + X(i') - X(f ) (261) 

Putting this into (258) we find 

v KC [n'](r, t) = J dt'd 3 r'W xc (n(r' - X{t'),t'), |r + X(i') - X(t)-r'| , t - t 1 ) 

= J dt'd 3 r'W xc (n(r'-Mt'),t')A(r-^(t))-(v'-X(t'))\,t-t') 

- v xc [n](r-X.(t),t). (262) 

Thus (258) satisfies (242). Hence the HPT and Galileian invariance requirements 
are satisfied. 
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Apart from its presumed nonlinear capabilities, (258), is also more general than 
the linear response work in preceding sections because it invokes the spatial nonlo- 
cality of the uniform gas xc kernel / xc , making it more comparable to the work of 
Dabrowski [102]. To compare it with the discussions above, we now make the local 
approximation 

t; r't') « S(v - r').f x h c OI >, q = 0, t - t') (263) 

If this local approximation is employed in Eqs. (256) and (257), the resulting 
t> xc [n](r, t) does not satisfy the HPT and the requirements of Galileian invariancc 
because, when linearized, it reduces to the Gross-Kohn form (191) and this is known 
[112] not to satisfy the HPT. However, combined with Eq. (258), the local approx- 
imation (263) leads to 

v xc [n](r,t) = J dt'w xc (n(R(t' \ r,t),t'),t - t') (264) 

where 

w xc (n, t) = / f*r(P, q = 0, r)dp. (265) 
Jo 

We now show that the functional (264), when linearized, gives precisely the modified 
linear-response xc kernel / xc of Eq. (248). To this end we consider small motions 
around a static equilibrium, in the sense that the displacement 

x(r,i) = R(i|r,i ) -r (266) 

of each fluid element from its initial (t — t n ) position r is small. Under these 
circumstances both the fluid displacement X and the current J in (259) are small 
(first-order) quantities. Thus we may write 

R(i'|r, t) = r + O(x) = R(t'|r, t ) + O(x) (267) 

and using this in (259) we find 



AR(t'|r, t) = J(R(t'|r, t ))/n (r) + 0(x 2 ) - ^Wl^ to) + 0(x 2 ). (268) 



Integrating both sides of (268) with respect to t' , starting from t' — i, we find 

R(t'|r, t) - R(i|r, t) = R{t'\r, t ) - R(t\r, t ) + 0(x 2 ), (269) 
so that, by (269) and (266) 

R(t' | r, t) = r + x(t') - x(i) + 0(x 2 ). (270) 
In the linear limit we can also integrate the linearized continuity equation 

^+V> o (r)u(r,t)]=0 (271) 

to give the density perturbation rii in terms of the fluid displacement x : 

ni(r,i) = -V • [no(r)x(r,t)] +0(x 2 ). (272) 

We can now use (270) and (272) to expand the density argument of w xc in the 
nonlinear functional (264): 

n(R(t'|r, t), t') = n(r + x(r, t') - x(r, t), t') + 0(x 2 ) 
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= n (r) + Vn (r) • (x(r, t') - x(r, i)) + m(r, t') + 0(x 2 ) 

= no(r) - Vn (r)-x(r, i) - n (r)V • x(r, t') + 0(x 2 ). (273) 

In deriving (273), we used (270) in the first step, standard linearization in the 
second step and (272) in the third step. Putting (273) into the proposed nonlinear 
xc potential (264) we find 

v xc (r,t)= w xc (n (r),t-t')dt' - Vn (r)-x(r,t) / -^(n (r), t - t')dt' 

(274) 

-Mr) I- x MO, t t')W ■ x(r, t>)dt> + 0(x 2 ) 
Integrating (256) with respect to time and using / x ° m (n) = dv^° m (n)/dn, we find 
that the first term in (274) is just the xc potential v^° m (n (r)) of the static, unper- 
turbed problem. The linear correction to this equilibrium value of the xc potential 
is then, by (274) and (256), 

/OO pOC 
fxc{n {r),t-t')dt')n 1B (r,t)+ / / xc (n (r), t - t')n 1A (r, t')dt'. 
-oo J —oo 

(275) 

where 

ni A (r,t) = -no(r)V-x(r,t), n 1B (r, t) = -x(r, t)-Vfi (r) (276) 
Fourier-transforming (275) and writing the terms in reverse order we find 

vi xc {r,u)) = f xc (n Q (r),u)n 1A (r,u)) + / xc (n (r),u = 0)n 1B (r,u)). (277) 

This is identical to the form (248). To summarize, we have proposed a rather bold 
Ansatz, Eq. (258), for the time-dependent xc potential v xc (r,i) of an arbitrary sys- 
tem which could be far from equilibrium. This Ansatz carries a nonlocal space and 
time dependence based on uniform-gas data, but accesses the actual system density 
n(r, t) in a simple local fashion. The assumption of local space dependence in the 
uniform gas yields a simpler form again, Eq. (264). It remains to be seen whether 
our relatively simple forms can cope with the gamut of nonlinear phenomena in 
systems far from equilibrium. As a first step it would be interesting to investigate 
second-and higher-order nonlinear susceptibilities described in section 5.2. Com- 
puter codes for investigating the fully nonlinear case may be adaptable from the 
work of Galdrikian et al. [123], who have investigated strongly driven quantum wells. 
Regardless of the applicability of our method to general nonlinear phenomena, the 
use of the trajectory function R(i' | r, t) in (258) and (264) guarantees two things: 
firstly, satisfaction of the generalized Galileian invariance condition (242) and hence 
of the Harmonic Potential Theorem (240), for motion of arbitrarily large amplitude; 
and secondly, for systems close to equilibrium the more local version (264) reduces 
to the linear time-delayed or frequency-dependent formalism previously proposed 
by Dobson [112]. 

6.2 Time-dependent optimized effective potential 

The approximate xc potentials described so far were derived from the homogeneous 
electron gas in one or another way. All of them have one deficiency in common: 
They contain spurious self-interaction contributions. It is known from static DFT 
that the removal of self-interaction is an important ingredient in the construction of 
good xc potentials. Various approaches to the construction of self-interaction-free 
functional are known in the static case [37, 124, 125, 126, 127, 128, 129, 130, 131, 
132, 133, 134, 135]. One of these is the so-called optimized potential method (OPM) 
[133, 134, 135]. This method takes as starting point a given expression for the total 
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energy E[<f>\ . . . cfy^] of an TV-electron system as a functional of a set of single-particle 
orbitals {</>.,■ (r)} (e. g. the Hartree-Fock total energy functional in the exchange-only 
case). Then, the variationally best local effective potential is determined such that, 
when inserted in a stationary single-particle Schrodinger equation, it yields the 
set of N eigenfunctions (corresponding to the N lowest eigenvalues) that minimize 
E[4>\ . . . </>jv]- In practice, the full OPM scheme is computationally quite involved 
since it requires the numerical solution of an integral equation for u xc (r). As a 
consequence, complete OPM calculations have been performed mainly for problems 
where the potential is a function of a single variable, e. g. for spherically symmetric 
atoms [134] -[140]. There exists, however, an approximate OPM scheme, recently 
proposed by Krieger, Li, and Iafrate (KLI) [141] - [149], which is numerically as 
easy to handle as the ordinary KS scheme. This simplified OPM has been applied 
very successfully to the calculation of atomic properties [6] . In many respects this 
method is currently the most accurate density-functional scheme. 

In the present section we shall describe the construction of a self-interaction- 
free xc potential which can be viewed as a time-dependent version of the optimized 
potential method (TDOPM). The approach leads to function of (r, t) rather 

than to v xc as an explicit functional of the density. In order to derive such a 
time-dependent generalization of the OPM we consider an TV-electron system at 
some finite time to which, for all times up until to, has been in the ground state 
associated with an external potential v (r) (e.g., a nuclear Coulomb potential). We 
assume that the corresponding stationary OPM problem has been solved for that 
system, i. e. a local effective potential and a set of N single-particle orbitals {4>j} 
(with energy eigenvalues Sj) minimizing a given energy functional E[<j>\ . . . 4>n] are 
assumed to be known. Again, at t = to an additional time-dependent potential 
Vi(r,t) is switched on. Our goal is to determine the time evolution of the system 
under the influence of the total external potential v(r, t) — v (r)+vi(r, t) from to up 
until an arbitrary later time t\. To construct an optimized local effective potential 
we start with the quantum mechanical action [150] 



JV „ t , , „ 2 



A[ipi...ip N ] = J2 J dtjd 3 rip*{r,t)(i 



j 

ft 



d_ 



<Pj(r,t) 



J 1 dtj d?r n(r, t)v(r, t) \ j ' dt j ' d\ J d 3 r> ^^'^ Ac [<pi • • . ^](278) 

written as a functional of N time-dependent single-particle orbitals {<fij{r, t)} where 

n(r,t) = J2j IVj( r )*)| 2 - I n a time-dependent exchange-only theory Ax C [<Pi ■ ■ ■ ¥n] 
— the xc part of the quantum mechanical action — would be replaced by the time- 
dependent Fock expression 

A„ - A - -if>„, ^Mll^MM (279) 

(<jj denotes the spin orientation of the jth orbital). We note that the integrand 
of (279) is a local expression with respect to the time-coordinate, i.e., all orbitals 
depend on the same time argument t. With approximate functionals of this type, 
the causality problem described in section 5.2 does not occur. The orbitals are 
solutions of the time-dependent Schrodinger equation 

.3 , , /V 2 



l dt 



<p j (r,t)= (-^-+v„(T,t)^ip j (T,t) , j = l,...,N , (280) 

with <fij(r, t) = (j>j(r) exp[— iej(t — to)] for — oo < t < to- The local effective potential 
is given by 

^(r,t)=«(r,i) + « H (r,t) +U J c DOPM (r,t) , (281) 
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where v H (r,t) = Jd 3 r'n(r', t)/\r— r'| denotes the time-dependent Hartree potential. 
The total potential v s (r, t) has to be determined in such a way that the {<Pj(r, t)}, 
resulting from Eq. (280), render the total action functional A[<pi ■ ■ ■ pn] stationary. 
Therefore, we have to solve the following variational problem: 

5A[<pi ...<p N ] ^ [ + °° ,., f, 3 , ( SA[tp! ...tp N ] 8p 3 (r',t') 5A[ip\ . . . ip N ] Sip* (r', t') ' 
" L L i 8 Vj {v',t') Sv s (r,t) + 



8v s (r,t) ^j-oo J 1 6i Pj( r '^) Sv *(r,t) 5<p*(r',t>) Sv s (r,t) J 

3 \ J / 

= . (282) 
We first compute the functional derivatives 5A/Spj and 5 A/ Sip*: defining 

1 SA xc {ipi . . .ip N ] 



= ^7T* ' ( 283 ) 



ip*(r,t) Sipj{v,t) 



we obtain 
5A[<pi . . . <p N ] 



~ l W~{~^2~ + W(r '' ^ + Wli(r '' ^ + Uxcj (r '' t,} ) 



(284) 

and an analogous expression for 5A/S<p* which, for all reasonable (i. e. real) func- 
tional A[p>\ . . . <P>n\, is the complex conjugate of (284). 0{x) denotes the usual step 
function (1 for x > 0, for x < 0). To arrive at Eq. (284) the first term of Eq. (278) 
has to be integrated by parts with respect to the time coordinate. We impose the 
usual boundary condition on ipj(r,t) at t — t\, i. e. 5ipj(r,ti) = 0, thus obtaining 
a zero boundary contribution. The other boundary contribution at t = — oo van- 
ishes, too, because the action functional (278), in order to be well-defined, is to be 
calculated by introducing the usual factor e vt in the integrand and taking lim^^ 0+ 
after the integration. Substituting Eq. (281) into (284) and making use of the fact 
that ip* solves the complex conjugate of the Schrodinger equation (280), we find 

A ll] {r ;^ ] - K T c D ° PM (r',0 - u xcj (r',t')] <p*(r',t') *( tl - t') . (285) 

In order to evaluate SA/Sv s from Eq. (282), we further need the functional deriva- 
tives S(pj/Sv s and Sip*/Sv s . The stationary OPM eigenfunctions {0j(r), j = 1, . . . , oo} 
form a complete orthonormal set, and so do the time-evolved states {ipj(r,t), j = 
1, . . . , oo} for any time t € [— oo, t{\ : and we denote this set by <f> t - Now consider $ 4 
as unperturbed states, remembering that at t = t\ the orbitals are held fixed with 
respect to variations in the total potential. We therefore start from t = ti, subject 
the system to an additional small perturbation 5v s (r, t) and let it evolve backward 
in time. The corresponding perturbed wave functions ip'j(r,t) are determined by 
the backward Schrodinger equation 

ij^.(r,t)= (- S ^-+v s (r,t)+Sv s (r,t) S J<p' j {r,t) , j = 1, ... ,N (286) 

with the initial condition <p'j(r,ti) — (pj(r,t\). This problem cannot be treated 
directly with time-dependent perturbation theory as described in standard text 
books because the unperturbed Hamiltonian is already time-dependent. Neverthe- 
less, Dirac's method of variation of constants can be applied in a straightforward 
manner. We expand, at each given t, the perturbed wave function ipj(r,t) in terms 
of the set $ t , 

oo 

^ j {T,t) = ^2c jk (t)( Pk (r,t) , (287) 
fc=i 
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and insert this expansion in (286), utilizing Eq. (280). The resulting equation 

oo oo 

i Cjfc(%fc(r, t)=^2 Cjk(t)Sv s (r, t)tp k (r, t) (288) 



i 

k=l fc=l 



is then multiplied by tp*(r,t) and integrated over all space; the orthonormality of 
$i yields 

1 00 f 

Cjl(t) - -J2c jk {t) d 3 r tf(r,t)Sv 8 (r,t)<p k (r,t) . (289) 
1 fe=i J 

We now make the usual ansatz for a perturbation expansion, 

c jk (t) = c<$(t)+c$(t) + ... (290) 
and collect corresponding orders on each side of Eq. (289). This yields 

cf(t) = 

= T^cg'W d*r y*{r,t)5v s {v,t)y k {v,t) (291) 

k=l 



Since, in our case, the wave function evolves backward from the fixed state <fj(r, t\) 
we find cf k \i) — &jk and c^(ti) — 0, leading to 

ft 



cf{t)= l J^dt'jd i ry* l {r,t l )8v s {r,t l )y j {v,t') . (292) 



It follows that the first-order correction to the wave function ipj(r,t) under the 
influence of Sv 8 (r, t) is given by 

S Vj {T, t) = jr, c? k \t)tp k (r, t) = i £ tdt' f dV^r', t')6v s (r>, t') Vj {v\ t> fc (r, t) . 
fc=i k=i Jt J 

(293) 

Therefore, the desired functional derivative is 

S sf(rt) = 1 £ ^ t)iPj (r ' * Vfe(r '' ^ 9{h t] 6{t ^ ■ (294) 
s ^ ' ' fe=i 

Once again, Sip* /5v s leads to the complex conjugate expression. We can now insert 
(285) and (294) in the variational equation (282), and the result is the TDOPM 
integral equation for the local exchange-correlation potential v KC (r,t): 

N rti r 

tJ2 / dt' d 3 r' [v^ OPM (r\t')-u xc] (r\t')}^(r,t) V *(r\t')K(r 7 tyX)+c.c. = 
j J— 00 J 

(295) 

The kernel 

00 

K(r, t, r', t') = <Pk(r, t)<PkV , *(* - *') ( 2 96) 
fc=i 

can be identified with the Green's function of the system, which satisfies the differ- 
ential equation 



7/2 



. I -(4- + ».(-',*') 



Jf(r, i, r', t') = -*<J(r - r')<*(* - i') (297) 
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with the initial condition K(r,t,r',t') = for t' > t. The TDOPM scheme is 
now complete: the integral equation (295) has to be solved for v xc (r,t) in combi- 
nation with the Schrddingcr equation (280) and the differential equation (297) for 
K(r, t, r', t'), both with the appropriate initial conditions. It is easy to show that 
in the time interval [— oo,ii] the exchange-correlation potential v xc (r,t) is only de- 
termined up to within a purely time-dependent function c(t) (as expected in view 
of the time-dependent HK theorem discussed in section 2). 

We now demonstrate that for t < t n or for a time-independent external potential 
(vi(r,t) = 0) the TDOPM reduces to the stationary OPM. For this purpose we 
rewrite Eq. (295) in the following way (using the fact that v xc is real) : 

JL ft! 

I 



£ dt' fd s r'[v^ OPM (r',t')-u xcj (r',t')] Vj (v,t)^(r\ + c.c. 

k^j 

*2^-(r,t)^(r,t) / d 3 r'(^(r',f')-< CJ (r',i'))^(r',f')^(r',f') . (298) 



In the static case, the orbitals {<Pj(r, t)} are replaced by {<j>j(r) exp[— iej(t — to)]}. It 
is reasonable to assume that the exchange-correlation functional A xc then becomes 

ft! 

A xc [ip!...ip N ]^ dt' E xc [<p 1 (t')...<p N (t')} , (299) 

J — oc 

where E xc [(j>i . . . 4>n] is the corresponding ground state exchange-correlation energy 
functional. Definition (283) then yields 



<c? IC M) = 



1 5E xc [4>i . . . cj) N ] 



(300) 

^(r)=^(r)e- iE i (t -'o) 



We assume that the value of E xc [cf>i ... 0^] remains unchanged if the arguments 
{(j)j(r)} are multiplied by phase factors e Mj '. If this is the case, we can use the 
identity 

(where ro is an arbitrary reference point) and write E xc in Eq. (299) as a functional 
of the combinations ipj(r 7 t)ip*(r' ,t). Then it is not difficult to show that u x ^ lc is 
independent of time and that the right-hand side of (298) is zero. We therefore 
obtain 

N ft! f 00 

]T / dt' dV K° c PM (r') - <f c (r')] U^ir') £ ^r^Oe"^-^-''^* t') 

k^j 

+ c.c. = . (302) 

Performing the integration over t' we find the stationary OPM integral equation 
[134] 

lim jr [ dV K° c PM (rO ~ < c f V)] (r)ff (r*) f] f^f^ + c.c. = . 

Mi 

(303) 

The derivation of Eq. (303) shows that in order to recover the static limit from the 
time-dependent formalism one had to extend the time integral in Eq. (278) to — oo; 
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a finite lower time boundary does not correctly account for memory effects in v xc 
and therefore results in an unphysical time dependence even in the static case. 

The numerical implementation of the full TDOPM is an extremely demanding 
task. It is therefore most desirable to obtain a simplified scheme. To this end we 
shall perform a transformation of Eq. (295) similar to the one proposed by KLI in 
the stationary case [146, 149]. This will lead to an alternative but still exact form 
of the TDOPM scheme which allows one to construct approximations of v xc (r,t) 
which are explicit functionals of the orbitals {<fj}, thereby avoiding the need to 
solve the integral equation. Following Refs. [146, 149], we define 



Pj(r,t) = 



Vi(r,*) 



df/dV [»r M (r',i')-^(r',/)]^(r',i')E^(M)^(r',0^ 

(304) 



k=l 
k^j 



and 



u xcj (*) = / d 3 r rij (r, t)u xcj (r, t) 



(305) 



where rij(r,t) = \<Pj(r, t)\ 2 . Eq. (298) can then be written as 



N 



N 



^2/fij{T,t)pj{r,t) + ex. = -i^Tnjfat) [ dt (u xcj (t) -u* xcj (t)) , 

(306) 



and it is easy to show that 



d 3 r nj(r,t)pj(r,t) = . 



(307) 



Evaluating tpj(r,t)[—i d/dt + V 2 /2 — v s (r,i)]<Pj(r,t)pj(r,t) one finds that pj(r,t) 
satisfies the following differential equation: 



1 d 

-V-(nj(r,t)Vpj-(r,t)) - in 3 (r,t) ^Pj(r.t) - i Jj(r, t) ■ Vpj(r, t) 
= -nj{r,t) [«™ OPM (r,t) - u xcj {r,t) - (v xcj {t) - u xc] (t) 



(308) 



with the current density Jj(r, t) = (2i) 1 (<p*(r, t)V<Pj(r, t) — <fj(r, i)Vtpj(r, t)) and 
v xcj (t) = jd 3 r nj(r,t)v™ OPM (r,t). Finally, operating with V 2 on Eq. (306) and 
using Eq. (308) we find 

..TDOPM / 



l (r,t) 



+ 



+ 



1 N 1 
— — ]T nj (r, t) - « cj (r, t) + < c . (r, t)) 

j 

1 N 

^)E^( r >*) 



4n(r 



V 2 n, (r, t) J Jt> (u XCJ (t>) - u* cj (t')) (309) 



where 

< C jM) = u xcj (r,t)+ 



rij(r,t) 



1 d 

-V-(pj (r, t) Vnj (r, t))+inj (r, t) — Pj (r, i)+iJj (r, t) -Vpj (r, t) 

(310) 



Eqs. (309) and (310) together with the differential equation (308) for pj(r,t) and 
the condition (307) (which can be used to fix the constant left undetermined by 
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Eq. (308)) represent an exact alternative formulation of the TDOPM scheme. The 
advantage of Eq. (309) lies in the fact that it is a very convenient starting point for 
constructing approximations of v xc (r, t) as explicit functional of the {<pj (r, t)}: it is 
only necessary to approximate Pj(r,t) in Eq. (310) by a suitably chosen functional 
of the orbitals. We can then readily solve Eq. (309) analytically for v™ OPM (r,i), 
as we shall show below. 

We expect an approximate potential v xc (r,t) defined in this way to be close to 
the exact v xc (r, t). This conjecture is based on the observation that the difference 
between v xc and v xc is entirely accounted for by the differences u' xc j — u XC j which 
are zero if averaged over the jth orbital, as will be demonstrated in the following. 
From Eq. (310) we obtain 



ljd 3 rV- { Pj (r,t)V nj (r,t))+iJ d 3 r 



n 3 ( r > *) g- t Pj ( r > *) + Jj ( r . t) ■ Vpj (r, t) 



311) 



u 



Using the divergence theorem, the first term on the right-hand side can be trans- 
formed into a surface integral which vanishes if the time-dependent orbitals decrease 
exponentially for r — > oo. The contribution to the second integral containing Jj-Vpj 
is then integrated by parts. The surface term vanishes due to the same argument 
as before, and the remaining term is transformed using the continuity equation for 
the jth orbital to replace —V • Jj(r, t) by drij(r, t)/dt. Hence we find 

Q n 

xcj (i) - u xc] (t)=ig- t J d 3 r nj (r, t)p 3 (r, t) = , (312) 

where the last equality follows from Eq. (307). 

The simplest approximation is obtained by replacing pj by its average value, 
i. e. by setting Pj(r, t) = 0. The resulting approximate potential will be termed the 
time-dependent KLI (TDKLI) potential. It is given by the equation 

1 N 1 

3 

1 N 

+ I^T7)E V2 ^M) J Jt'(u xcj (t')-K cj (t')) -(313) 



TDKLI 



vl^ l (t)-\(u^(t)+ul c M 



This equation is still an integral equation for u™ KLI . It can, however, be solved 
semi-analytically [145]: Multiplying Eq. (313) by rife(r, t) and integrating over all 
space yields 



JV 

-TDKLI 

J xck 



(t) = w xck (t)+Y J M kj (t)vl^\t) , (314) 



where we have defined 



1 N 1 

w xc (r,t) = -t-tyE 71 ^ 1- '^ o ( M ^-(r,t) +u* cj (r,t)) 



n(r,t) 

(*)) 



3 
N 



n(r,t) 

+ ^J) E v2 "i ( r ' *) / J** 7 (*') - (O) (315) 
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N 



4G 



and 

Solving Eq. (314) for «™ KLI (i) requires inversion of the TV x N matrix 

A kj {t) = S kj - M kj (t) (317) 

and leads to 



N 

v™ KL1 (t) = Y / (A- 1 (t)) jk W xck (t) . (318) 



When Eq. (318) is substituted into Eq. (313), one obtains v™ KLI (r,t) as an ex- 
plicit functional of the orbitals {(fj(r,t)}. As the exact v xc (r,t) which follows from 
Eq. (295), v™ KLI (r,i) is determined by Eq. (313) only up to within a purely time- 
dependent function c(t). 

The last term of Eqs. (313) and (315) vanishes identically for a large class of 
exchange-correlation functionals A xc - This class includes all functional depending 
on {fj} only through the combinations ipj(r, t)(pj(r', t) (such as the time-dependent 
Fock functional, Eq. (279)). 

One readily verifies that both the full TDOPM potential and the TDKLI ap- 
proximation of it satisfy the generalized translational invariance condition (242) 
(and hence the harmonic potential theorem) provided that 

Acb'i • ■ • <p' N ] = Acbi • • ■ <Pn] (319) 
is satisfied with ip'j being the orbitals in the accelerated frame: 

$(r,t) =exp(-iS(t) + iXT)0 j (r-X(t)). (320) 

The TDHF functional (279) is easily seen to satisfy the constraint (319). Equation 
(319) will be a strong guideline in the proper construction of approximate correlation 
functionals A c [(fii . . . <pn\- Equation (313) combined with the Schrddinger equation 
(280) represents a time-dependent scheme which is numerically much less involved 
than, e.g., the time-dependent Hartree-Fock method. Numerical results obtained 
with this scheme for atoms in strong laser pulses will be described in section 8. 

To conclude this section we construct in the following an approximation of the 
xc kernel f xc on the basis of the TDOPM. A calculation analogous to Eqs. (138) 
- (152) shows [151, 152] that within TDOPM the central equation (152) holds for 
the quantity /™ OPM defined by the integral equation 



^(r)0;(r> fc (r')#(r)e-*< £ ' 



x^fjfl^iy^r'^-g^iy,^',^) c.c. 



(321) 



where 



i i,,TDOPM/ v _-\ 

flO-)( v r r' f') = 

' ly "-' - ' 2^(r',t') 5 Vj (r>,t>) 



(322) 

Equation (321) has the same algebraic structure as the TDOPM integral equa- 
tion (295) with the time-dependent orbitals tpj(r,t) replaced by e _l£j *(/>j an d with 

« x T c DOPM (r',t') and u xci (v',t') replaced by / x T c DOPM (y,r,r',t') and <${y,T,r> ,t), 
respectively. A simple analytical approximation to v J c DOPM (r, t) is given by 

<P r ° X (M) = E ' 2n(r j)' 2 KcJ-fr*) +<c 3 (r,*)) • (323) 
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Applying this approximation to (321), i.e., setting 



/ x a r ox (y, r, r', t') = J2 J fSf (gg (y, r, r', f ) + c.c.) (324) 

j 

and using the explicit analytical form (323) to evaluate (322) one arrives in the 
time-dependent x-only limit (279) at the compact expression 

,/ xc InoKr.r.w)- | r _ r ,| no(r) no(r ,) • 

In general, the Fourier transform of the xc kernel defined by Eq. (321) is frequency 
dependent (even in the TD x-only case), a feature which is not accounted for by the 
present approximation (325). However, for the special case of a two-electron system 
treated within TD x-only theory, Eqs. (323) and (325) are the exact solutions of the 
respective integral equations. 



7 Applications within the perturbative regime 

7.1 Photoresponse of finite and infinite Systems 

To date, most applications of TDDFT have been in the linear response regime. 
Calculations of the photoresponse from Eqs. (154) and (155) are, by now, a mature 
subject. The literature on such calculations is enormous and a whole volume [153] 
has been devoted to the subject. In this section we shall restrict ourselves to the 
basic ideas rather than describing the applicational details. 

The TDKS formalism has been employed to calculate the photoresponse of atoms 
[10, 12, 13, 14, 154, 155], molecules [156, 157] and clusters [158, 159, 160, 161, 162, 
163, 164, 165, 166, 167, 168] metallic surfaces [169, 170, 171, 172, 173, 173, 174, 175] 
and semiconductor heterostructures [72, 176, 177, 178, 179] bulk semiconductors 
[180] and bulk metals [181, 182, 183, 184] 

For simplicity, we consider sufficiently low radiation frequencies, such that the 
electric field can be assumed to be constant across the atom or molecule. For atoms 
this is the case for photon energies below 3 keV. The external potential associated 
with a monochromatic electric field is then given by 

t>i(r,t) = Ezcoswt. (326) 

The induced density change n(r, t) — n (r) (143) can be characterized by a series of 
multipole moments. The induced dipole moment 

p(t) = - J d 3 rz(n(r, t) - n (r)) , (327) 

can be expanded as [185] 

p{t) = a(u>) ■ Ecoswt + i/3(0):EE+ i/3(2w):EEcos 2wi (328) 

+ |7(w):EEEcoswt + i7(3w):EEEcos3wt, (329) 

where the notation is meant to indicate the tensorial character of the quantities a, (3 
and 7. The first coefficient, a, is termed the dipole polarizability; (3 and 7 denote 
the second-and third-order dipole hyperpolarizabilities. For spherically symmetric 
ground states (3 is zero. 
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Figure 3: Total photoabsorption cross section of the Xe atom versus photon energy 
in the vicinity of the 4d threshold. Solid line: self-consistent time-dependent KS 
calculation from [10]; crosses: experimental data from [186]. 

The dipole polarizability is related to the frequency-dependent linear density 
response n\ (r, u) via 



Zangwill and Soven [10] have calculated the photoabsorption spectrum of rare-gas 
atoms from the frequency-dependent KS equations (156) - (157) within the ALDA. 
As an example for the quality of the results we show, in Fig. 3, the photoabsorption 
cross section of Xenon just above the 4d threshold. The agreement with exper- 
iment is remarkably good. Results of similar quality have been achieved for the 
photoresponse of small molecules [156, 157]. 

It should be mentioned that the thresholds characterizing the onset of continuous 
absorption from the various occupied atomic shells are not well reproduced in the 
calculations of Zangwill and Soven. The calculated absorption edges are typically 
several eV below the observed thresholds. While, in principle, TDDFT should yield 
the correct thresholds, it appears that simple approximations such as the ALDA 
are not sufficient in this respect. 

As a point of practical interest we mention that the KS response function is 
usually not calculated directly from the KS orbitals as in Eq. (157). Instead, one 
rewrites the response function in terms of the KS Green's function. The latter is 
then calculated numerically from the corresponding equation of motion, usually by 
multipolc expansion [10, 187]. 

The linear photoresponse of metal clusters was successfully calculated for spher- 
ical [158, 160, 159, 163] as well as for spheroidal clusters [164] within the jellium 
model [188] using the LDA. The results are improved considerably by the use of 
self-interaction corrected functional. In the context of response calculations, self- 
interaction effects occur at three different levels: First of all, the static KS or- 
bitals, which enter the response function, have a self-interaction error if calculated 
within LDA. This is because the LDA xc potential of finite systems shows an ex- 
ponential rather than the correct — 1/r behaviour in the asymptotic region. As a 




(330) 



and the photoabsorption cross section is given by 

. . 47TO> , 
<T(L0) = Xsa[LLl) . 



(331) 
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consequence, the valence electrons of finite systems are too weakly bound and the 
effective (ground-state) potential docs not support high-lying unoccupied states. 
Apart from the response function \s , the xc kernel / xc [no] no matter which approx- 
imation is used for it, also has a self-interaction error. This is because / xc [^o] is 
evaluated at the unperturbed ground-state density n (r), and this density exhibits 
self-interaction errors if the KS orbitals were calculated in LDA. Finally the ALDA 
form of / xc itself carries another self-interaction error. 

To improve upon these defects, one has to go beyond the LDA: The (modified) 
weighted density approximation [189] retains the correct asymptotic behaviour of v xc 
and improves the response properties of metal clusters [162, 165]. A different route 
to improvement provides the self-interaction correction (SIC) of Pcrdcw Zunger [37], 
where the spurious self-interaction of the LDA is compensated by additional terms 
in the ground-state potential [166] and in the effective perturbing potential as well 
[167] (Full-SIC). 

In most theoretical work on the response of metallic surfaces the ionic potential 
is replaced by the potential due to a uniform positive charge background in a half 
space, say z > 0. This is the so-called jellium model for metallic surfaces. In this 
model are two intrinsic microscopic length scales, the inverse Fermi wave-number, 
kp 1 , and the Thomas-Fermi screening length (« surface thickness), k^ F . Both 
lengths are typically of the order a w 10~ 8 cm. In most applications the perturbing 
electric potential V\ and the perturbing vector potential A x vary on a length scale 
i which satisfies i 3> a. Examples arc the scalar potential t>i(r) due to an external 
charge at a distance z 3> a, or the vector potential Ai(r, t), associated with a 
light wave of wavelength A 3> a. The corresponding linear responses n\ and ji 
vary on the scale of I in the x — y plane but, because of the abrupt drop of the 
unperturbed density at the surface (on the scale of a), they vary on the short scale 
a in the z-direction. Formal arguments due to Feibelman [190] have shown that, to 
leading order in a/£, the effect of the surface on the electromagnetic fields far from 
the surface (\z\ >• a) is entirely characterized by two complex frequency-dependent 
lengths, d||(w) and d±(u). 

Actually, for the jellium model du (u>) = 0. This result has been obtained in 
the random phase approximation (RPA) in Ref. [190]. It is easily established as a 
rigorous many-body result for the jellium model [191]. To define d±(ui) we Fourier 
analyze all physical quantities parallel to the surface, in the x— y plane. For example, 
a Fourier component of the induced charge density becomes 

ni(r,w) =ni(z,cj)e 4q ii' r , (332) 

where = (q x , q y , 0) (and |q|| | a <C 1). Then d±(iv) is given by 

d±{uj) = J - T - -, (333) 

J dzni(z,LU) 

i.e., it is the (complex) center of mass of the induced surface charge. d±(u) is the 
generalization of the static image plane introduced by Lang and Kohn [192]. 

This d±(ui) is then the subject of quantitative calculations. They require the den- 
sity response n\(z,u)), to a uniform external electric field perpendicular to the sur- 
face. The calculation was first carried out in the RPA equivalent to time-dependent 
Hartree theory, in which the xc kernel / xc is neglected. These calculations led to very 
interesting results not present in classical Maxwell theory, such as the surface photo 
effect and surface plasmons. Plasmons are high-frequency charge-density oscilla- 
tions of the electron gas. In a bulk material the long-wavelength plasma frequency 
is ujp — (Anne 2 /m) 1 / 2 in gaussian cgs units. Plasmons occur in the ultraviolet 
frequency region for metals, but the artificial electron gas in semiconductor quan- 
tum wells often has a plasma frequency in the infrared. The confinement of the 
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electron gas at metal edges introduces a range of new plasmon modes at frequen- 
cies other than uip, and these could potentially yield information about inhomoge- 
neous xc effects. Information about plasmons on films, surfaces and semiconductor 
wells is most easily available experimentally for small values of the surface-directed 
wavenumber q\\, and unfortunately in this region there are theorems prescribing 
the plasmon frequencies, regardless of the effects of exchange and correlation. A 
summary of these "no-go" theorems is given in [113] and further review is given in 
[119] for the case of semiconductor quantum wells. 

Early theoretical studies of plasmas used hydrodynamics [193]. These treatments 
were able to predict the main new feature of the plasmon spectrum at a metal surface 
due to the strong surface inhomogencity of the electron gas, namely the surface 
plasmon. Its frequency approaches ujp/\/2 as the surface directed wavenumber q\\ 
approaches zero, and this is correctly predicted in hydrodynamic and microscopic 
theories. This result is independent both of the precise electron density edge profile 
and of the type of xc kernel used, if any [194, 195]. Thus, although the surface 
plasmon is often the strongest feature in electron energy loss measurements on 
thin metal films, [196] it is hard to obtain any information from it about dynamic 
exchange and correlation. To see such effects one needs to measure with great 
accuracy the dispersion of surface plasmons. Only in the last few years has it 
been possible even to confirm experimentally a result first predicted by Feibclman 
[190] on the basis of selfconsistent RPA calculations, namely that the dispersion 
of the surface plasmon on a charge-neutral metal surface is initially negative. This 
result follows basically from the very "soft" or weakly bound nature of the electron 
gas at a neutral jellium surface, allowing electrons to spill out substantially into the 
vacuum. For a review of some experimental and theoretical aspects see [197]. While 
the value of this negative dispersion coefficient does depend to a degree on the xc 
kernel i xc introduced earlier, it remains to be seen whether experiments on metal 
films and surfaces can measure this quantity to a useful accuracy. On the theory 
side, an important observation by Liebsch [198] is that the KS orbitals used to 
construct the dynamic response must come from a static calculation using a model 
of exchange and correlation that is consistent with the dynamic xc kernel used in 
the plasmon calculation. For example, LDA calculation followed by RPA screening 
(with f xc = 0) is not consistent and causes false shifts in predicted surface plasmon 
frequencies. 

The weak binding and wide inhomogeneous density layer at the edge of a neutral 
metal surface leads to a "multipole" surface plasmon mode in addition to the usual 
"monopolc" surface plasmon [173, 172, 197]. This mode is in principle sensitive to 
f xc even at q\\ = 0. Gics and Gerhardts [173] and Dobson and Harris [199] investi- 
gated this mode both in the ALDA and the frequency dependent parametrization 
(206) - (210). It was found that, for an aluminium surface, the inclusion of the fre- 
quency dependence of / xc . has only a 3% effect on the multipole plasmon frequency, 
but a 20% effect on the damping of the mode. It seems likely that the frequency 
dependence of / xc will have a much larger effect on this mode for a low-electron 
density metal such as Rb, and this may be worth pursuing. 

In general, low-dimensional, low-density systems offer the best prospects for 
strong effects of xc phenomena on plasmon frequencies. A case in point is a pair of 
parallel quasi-two-dimensional electron layers in a semiconductor double-quantum 
well experiment. Interesting effects are predicted for this case [200]. 

Another way of probing dynamic xc effects experimentally is by inelastic X-ray 
scattering from bulk metals [201, 202, 203]. In this way, the so-called dynamical 
structure factor S(q,cj) can be measured which is proportional to the imaginary 
part of the full response function in reciprocal space. With this information at hand 
and with a first-principles calculation of the non-interacting response function, the 
connection (159) between / xc and the response functions can be used to deduce 
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information about / xc [204]. 

All applications quoted so far were for the linear response. Very few investi- 
gations have dealt with the higher-order response described in section 5.2. The 
frequency-dependent third-order hyperpolarizabilities of rare-gas atoms were calcu- 
lated by Senatore and Subbaswamy [86] within the ALDA; the calculated values 
turned out to bee too large by a factor of two, further indicating the need for self- 
interaction corrected functionals in the calculation of response properties. The effect 
of adsorbates on second-harmonic generation at simple metal surfaces was invested 
by Kuchler and Rebentrost [205, 206]. Most recently, the second-order harmonic 
generation in bulk insulators was calculated within the ALDA [207] . 



7.2 Calculation of excitation energies 

The traditional density-functional formalism of Hohenberg, Kohn and Sham [1, 2] 
is a powerful tool in predicting ground-state properties of many-electron systems 
[3, 4, 5]. The description of excited-state properties within density-functional the- 
ory, however, is notoriously difficult. One might be tempted to interpret the Kohn- 
Sham single-particle energy differences ojjk ■= £j — e k as excitation energies. This 
interpretation, however, has no rigorous basis and in practice the Kohn-Sham or- 
bital energy differences u)jk deviate by 10-50% from the true excitation energies 
Cl m := E m — Eq. Several extensions of ground-state DFT have been devised to 
tackle excited states. They are based either on the Rayleigh-Ritz principle for the 
lowest eigenstate of each symmetry class [208, 209, 210] or on a variational principle 
for ensembles [211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222]. A fun- 
damental difficulty is that the xc energy functionals appearing in these approaches 
depend on the symmetry labels of the state considered or on the particular ensem- 
ble, respectively. Until today very little is known on how these excited-state xc 
functionals differ from the ordinary ground-state xc energy. 

In this section we are going to develop a different approach to the calculation 
of excitation energies which is based on TDDFT [69, 84, 152]. Similar ideas were 
recently proposed by Casida [223] on the basis of the one-particle density matrix. 
To extract excitation energies from TDDFT we exploit the fact that the frequency- 
dependent linear density response of a finite system has discrete poles at the exci- 
tation energies of the unperturbed system. The idea is to use the formally exact 
representation (156) of the linear density response ni(r, w), to calculate the shift 
of the Kohn-Sham orbital energy differences ujjk (which are the poles of the Kohn- 
Sham response function) towards the true excitation energies fl m in a systematic 
fashion. 

The spin-dependent generalization [59] of TDDFT described in section 4.1 leads 
to the following analogue of Eq. (156) for the linear density response of electrons 
with spin a: 

n la (r,w) = / rf3 M s ^(rj;wK(y,w) (334) 

v J 

+ ^2 j d3 y J d3 y'xsa V (r,y;oj) ( j y ^ y /| +/xc^'(y,y';^ n lv ,{y',u). 

Here the spin-dependent exchange-correlation kernel is given by the Fourier trans- 
form of 

f (r- + r-> y-M — ^xcaK,nJ(r,t) 

,/xc GO 1 V 1 j 1 i ^ ) • 



(335) 



<W(r',t') 

with respect to (t — t'). Note that the spin-dependent response- function of nonin 
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teracting particles 



Xs*a'(r,r ;u>) = <W }_^(fk* - fja) ( _ g , (336) 

. , W l fc jtx fc fccr ) T *'/ 

is diagonal in the spin variable and exhibits poles at frequencies ujjka = tja — ^ka 
corresponding to single-particle excitations within the same spin space. In order to 
calculate the shifts towards the true excitation energies f2 of the interacting system, 
we rewrite Eq. (334) as 



1 



y-y 



]T J d 3 y' (^J(r-y')-X; / ^Xsa,(r,y; 

I/' \ f 

+/xci/i/'(y,y';w))^ni^(y',a;) = ^ J d 3 y Xsa»(r,y;u)v lu (y,uj) . (337) 

Since, in general, the true excitation energies i7 are not identical with the Kohn- 
Sham excitation energies ujjka, the right-hand side of Eq. (337) remains finite for 
u> — ► Cl. In contrast, the exact spin-density response n 1(T , has poles at the true 
excitation energies u> = il. Hence the integral operator acting on n la on the left- 
hand side of Eq. (337) cannot be invertible for u — > O. If it were invertible one 
could act with the inverse operator on both sides of Eq. (337) leading to a finite 
result for to — > il on the right-hand side in contradiction to the fact that ni a , on 
the left-hand side, has a pole at w = O. 

The true excitation energies f2 can therefore be characterized as those frequencies 
where the eigenvalues of the integral operator acting on the spin-density vector in 
Eq. (337) vanish or, if the integration over the delta- function is performed, where 
the eigenvalues A(w) of 



E 



Zj d3 y' E / d 3 yxs.Ar,y;u) (^^ + /xc^'(y, y';^)) My» = A(w)CM338) 



satisfy 

A(Q) = 1 . (339) 

This condition rigorously determines the true excitation spectrum of the interacting 
system considered. 

To simplify the notation, we now introduce double indices q = (j, k) so that 
ui qa = ej a —eka denotes the excitation energy of the single-particle transition (ka — ► 
ja). Moreover, we define 

$„(r) := <j>k*(r)*<f> Ja {r), (340) 
a q<7 := fka - fja (341) 

and set 

CM :=E/ d VE / <* 3 2/<W<My)* (^^ + / X c^(y,y';^))c,'(y>). 

(342) 

With these definitions, Eq. (338) takes the form 

E JX^+L ^H = A(«)Wr,«) • (343) 
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Solving this equation for Cr(r, w) and reinserting the result on the right-hand side 
of Eq. (342) leads to 



a' q 

with the matrix elements 



(344) 



M ?ff9 v'M =tt q >*> Jd 3 r j dV $* CT (r) + / XCCT(7 ,(r, r'; w)^ $ g v( r ')- 

(345) 

Note that the summation in Eq. (344) extends over all single-particle transitions q'a' 
between occupied and unoccupied Kohn-Sham orbitals, including the continuum 
states. Up to this point, no approximations have been made. In order to actually 
calculate A(w), the eigenvalue problem (344) has to be truncated in one way or 
another. One possibility is to expand all quantities in Eq. (344) about one particular 
KS-orbital energy difference u pT : 



AM 



A(iOp T ) 



dui 



(U> - LU pT ) + ... 



+ B(oj pT ) + ... 



(346) 
(347) 



The matrix elements with (uj pT ^ uj q ><j') can be written as 



M qaq , a ,(uj) 



uj - uj q , a , + irj Lu pa - Lu q > a < + irj du: 
whereas if (ui pr — w q i a i), 



M q<Tq > a <(uj) 



uj - Lu ql<J > + irj 



(W-Wpr) + ... (348) 



M qaqlcr >(uu) _ M qaqlcrl (uj pT ) dM qaqlcrl (u) 



u) — u) a 



irj u) — iv pT + irj 



+ 



dui 



+ ... 



(349) 



Inserting Eqs. (346) - (349) in Eq. (344) the coefficients A and B are readily iden- 
tified. If the pole Lo pT is non-degenerate, one finds: 



M^pr) = Mp TpT (ujp T ) 



and 



B(lu pt ) = dMpTpT 



du> 



I 



q' a'^pr 



M, 



pr q' a' 



(ujp T )M q i a i pT {uJp T ) 



LU pT - U) q 'cr< + 17] 

The corresponding eigenvector (in lowest order) is given by 

I 



qa 



A(u p 



-M, 



qa pr {^pr^pr 



(350) 



(351) 



(352) 



with (pr) fixed. The number £ pr is free and can be chosen to properly normalize 
the vector £. 

If the pole ujp T is p-fold degenerate, 

w pm = LJ P2r 2 = ■ ■ ■ = u PpTp = uj , (353) 

the lowest-order coefficient A in Eq. (347) is determined by the following matrix 
equation 



M PiTi PkTk M4;l = M^pfn , i = 1 • • • P ■ 



(354) 



fe=i 
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In general, one obtains p different eigenvalues A\ . . . A p . Then the remaining com- 
ponents of the corresponding eigenvectors can be calculated from 

^ = A^Joj P^P^M^l, (355) 

once the eigenvalue problem (354) has been solved. Assuming that the true excita- 
tion energy Q is not too far away from luq it will be sufficient to consider only the 
lowest-order terms of the above Laurent expansions. In particular, we set 

A„M « . (356) 

The condition (339) and its complex conjugate, A*(fi) = 1, then lead to 

Sl n =wo + ZtA n (u ) (357) 

This is the central result of our analysis. Eq. (357) shows that a single KS pole 
can lead to several many-body excitation energies. The corresponding oscillator 
strengths can be obtained [152] from the eigenvectors and the KS oscillator 
strengths. 

In the following, we exclusively consider closed-shell systems. For these systems, 
the Kohn-Sham orbital eigenvalues are degenerate with respect to the spin variable, 
which implies a lack of spin-multiplet structure. In what follows, we demonstrate 
how this is restored by the lowest-order corrections (357). Assuming that there 
are no further degeneracies besides the spin degeneracy, Eq. (354) reduces to the 
following (2 x 2) eigenvalue problem: 

53 M papa ,{yj )£ pa ,{uj ) = A£, ptT {uj Q ) . (358) 

For spin-saturated systems, = M P i P i and M p ^ p ^ = M p ^, so that the eigen- 

values of Eq. (358) are given by 

Ai, 2 = M pM ± M pVpi . (359) 

By Eq. (357), the resulting excitation energies are: 

fii - wb + R{Afp T p T +Mp Tpi } (360) 
fi 2 = wb + R{M pTpT -M p]pi } . (361) 

Inserting the explicit form of the matrix elements (345) one finds 

fii - u; a + 2VlJ d 3 r J d 3 r'$;(r) + / xc (r, r> )) $ p (r') (362) 

fi 2 = w + 23fj J <?r j d 3 r'^;{r)^ G xc (r,r';uj a )%(r') (363) 

where, for simplicity, we have dropped the spin-index of & pcr . 4 Obviously, the xc 
kernel appearing in Eq. (362), 

/ xc (r,r» = ± ]T jw(r,r';u;) (364) 

<7,<T' = ±1 



4 This is possible only if the unperturbed KS ground-state determinant is spin-saturated since, 
in this case, <f>jf (r) = <fijl( r ) f° r ^31 j. 
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is identical with the one already defined in section 5.1. On the other hand, Eq. (363) 
exhibits the kernel 



This quantity gives rise to exchange and correlation effects in the Kohn-Sham 
equation for the linear response of the frequency-dependent magnetization density 
m(r,u)) [59]. The fact that the magnetization density response naturally involves 
spin-flip processes, suggests that S^2 represents the spin triplet excitation energies of 
many-electron systems with spin-saturated ground states. The corresponding spin 
singlet excitation energies, on the other hand, are given by f2i . This assignment 
will be given further evidence by the numerical results presented at the end of this 
section. 

Apart from the truncation of the Laurent series, two further approximations are 
necessary: 

(i) The frequency-dependent xc kernels / xc and G xc have to be approximated. 

(ii) The static Kohn-Sham orbitals entering Eqs. (362) and (363) (cf. Eq. (340)) 
have to be calculated with an approximate (static) potential i> x ' at . 

As an application of the method, we consider the lowest excitation energies of the 
alkaline earth elements and the zinc series. Here, in addition to the degeneracy with 
respect to the spin index, the s — > p transitions under consideration are threefold 
degenerate in the magnetic quantum number m of the "final" state. Hence, we 
have six degenerate poles and Eq. (354) is a (6 x 6) eigenvalue problem. In our 
case, however, the matrix M PiTiPkTk in Eq. (354) consists of (three) identical (2 x 2) 
blocks, leading only to two distinct corrections, independent of m, as it should be. 

Tables 1-3 show the results of calculations based on Eqs. (362) and (363). The 
calculation of Table 1 employs the ordinary local density approximation (LDA) for 
u x * at and the adiabatic LDA (188) for / xc (both using the parametrization of Vosko, 
Wilk and Nusair [90]). In this limit, the kernel G xc is approximated by [103] 



The exchange-correlation contribution to the so called "spin-stiffness coefficient" 
a xc is also approximated within the LDA of [90] . 

The calculation of Table 2 uses the x-only optimized effective potential (OPM) 
for i> x * at in the approximation of Krieger, Li and Iafrate (KLI) [224] and for / xc the 
TDOPM kernel (325) derived in section 6.2. Concerning the singlet spectrum, the 
OPM values are clearly superior to the LDA results and are also better than the 
usual Ascf values. The unoccupied orbitals and their energy eigenvalues are very 
sensitive to the behavior of the potential far from the nucleus. Thus one major 
reason for the superiority of the optimized effective potential is the fact that it is 
self-interaction free and therefore has the correct — 1/r tail (while the LDA potential 
falls off exponentially) . An important point to note is that the optimized effective 
potential decreases correctly for all orbitals. For this reason, the x-only optimized 
effective potential is also superior to the Hartree-Fock (HF) potential which is self- 
interaction free only for the occupied orbitals but not for the unoccupied ones. As 
a consequence, HF orbital-energy differences are typically too large. However, in 
spite of the fact that the OPM provides self-interaction free orbitals, it reproduces 
the triplet spectrum less accurately. This hinges on the approximation the xc kernel 
is based on. Substituting the TD-Fock expression (279) for the xc action functional 
defined in (278) leads to a xc kernel diagonal in spin space, because the correlation 




(365) 



.A LD A [n](r)r ^ )= , (r _ r0 _i_ 



a xc (n(r)) . 



(366) 
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between antiparallel spins is neglected. Accordingly, from Eqs. (364) and (365) we 
have, within the x-only TDOPM 

G x T c D ° PM W(r,r'^) = -l/™ OPM [n](r,r'; W ) . (367) 

This should be cured by adding appropriate correlation terms to the xc part 
of the action functional, which is further backed by the observation that when 
combining the advantage of approximating w x * at by the optimized effective potential, 
together with a local density prescription of exchange and correlation in the xc 
kernels / xc and G xc , both singlet and the triplet spectrum are reproduced well by 
Eqs. (362) and (363), as can be seen from Table 3. 

In spite of the fact that we focused our attention to the situation of closed shells, 
and spin-multiplets, the method is also capable of dealing with open-shell systems 
and spatial multiplets. More details can be found in [152]. 

We emphasize that the calculation of excitation energies from Eqs. (362) and 
(363) involves only known ground-state quantities, i.e., the ordinary static Kohn- 
Sham orbitals and the corresponding Kohn-Sham eigenvalues. Thus the scheme 
described here requires only one selfconsistent Kohn-Sham calculation, whereas the 
so-called Ascf procedure involves linear combinations of two or more selfconsistent 
total energies [209]. So far, the best results are obtained with the optimized effec- 
tive potential for u x * at in the KLI x-only approximation. Further improvement is 
expected from the inclusion of correlation terms [6, 225] in the OPM. 

7.3 Van der Waals interactions 

While TDDFT has its main applications in time-dependent phenomena, and in the 
calculation of excitation or promotion energies, certain aspects of groundstate energy 
calculations are also assisted by TDDFT. This development principally concerns 
the van der Waals (vdW) or dispersion-force component of the groundstate energy. 
The usual groundstate LDA and its various gradient extensions [227] do not give 
an adequate description of vdW forces [228] , presumably because these forces arise 
(in one picture at least: see below) from the correlations between dynamic electron 
density fluctuations in widely separated positions. This makes the usual local or 
near-local approximations invalid. The approach to be introduced here facilitates 
the derivation of van der Waals functionals via a frequency integration over dynamic 
susceptibilities. 

(ljvdW interactions for widely-separated fragments: Perhaps the most familiar 
example of a dispersion interaction is the attractive mutual energy of a pair of neu- 
tral spherical atoms separated by a large distance R, an interaction which forms the 
tail of the well-known Lennard- Jones potential. To lowest order this interaction en- 
ergy falls off [229] as R~ 6 . This form of dispersion interaction is readily derived for a 
general pair of non-overlapping electronic systems by regarding the electrons on the 
first system as distinguishable from those on the second system. One then obtains 
the R -6 dispersion energy (in addition to some "polarization " terms relating to 
any static electric moments [230]) by performing second-order Rayleigh-Schrodingcr 
perturbation theory, treating the Coulomb interaction between the two groups of 
electrons as the perturbation Hamiltonian . (For very large separations R the retar- 
dation of the electromagnetic interactions between the systems cannot be ignored. 
In this regime the R -6 law just quoted is replaced [229] by R~ 7 . This retarded form 
takes over whenever R» c/uj, where w is a characteristic response or fluctuation 
frequency of the electronic systems. We will consider only the non-retarded case 
here) . 

From the work of Casimir, Lifshitz, London and many others [229] we know 
that the perturbation expression for the dispersion interaction between separated 
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Atom 


State 


^cxp 


1 ISA 


^(A SC f) 


1 Cl A 


Be 


1 r> 
Pi 


0.388 


0.399 


0.331 


0.257 




Po 


0.200 










3 Pi 


0.200 


0.192 


0.181 


0.257 




3 P 2 


0.200 








A /Try 

IVlg 


1 P, 
"1 


u.oiy 


U.OOl 


n 9QQ 






3 Po 


0.199 










3 P1 


0.199 


0.209 


0.206 


0.249 




3 P2 


0.200 








Ca 


'Pi 


0.216 


0.263 


0.211 


0.176 




3 Po 


0.138 










3 Pl 


0.139 


0.145 


0.144 


0.176 




3 P2 


0.140 








Zn 


l Pl 


0.426 


0.477 


0.403 


0.352 




3 Po 


0.294 










3 Pl 


0.296 


0.314 


0.316 


0.352 




3 P2 


0.300 








Sr 


'Pi 


0.198 


0.241 


0.193 


0.163 




3 Po 


0.130 










3 Pl 


0.132 


0.136 


0.135 


0.163 




3 P2 


0.136 








Cd 


'Pi 


0.398 


0.427 


0.346 


0.303 




3 Po 


0.274 










3 Pl 


0.279 


0.269 


0.272 


0.303 




3 P2 


0.290 









Table 1: The lowest S-^P excitation energies of various atoms. The experimental 
values (Erst column) [226] are compared with results calculated from Eq. (362) for 
the singlet and from Eq. (363) for the triplet (second column) and with ordinary 
Ascf values (third column). The LDA was employed for v xc and the ALDA for the 
xc kernels. The corresponding Kohn-Sham orbital-energy differences are shown 
in the last column (All values in rydbergs). 
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Atom 


State 


^cxp 


n OPM 


^(Ascf) 


, 1VI 

uj 


Be 


1 r> 
Pi 


0.388 


0.392 


0.331 


0.259 




Po 


0.200 










3 Pi 


0.200 


0.138 


0.181 


0.259 




3 P 2 


0.200 








A /Try 

IVlg 


1 P, 
"1 


u.oiy 




n 9QQ 






3 Po 


0.199 










3 P1 


0.199 


0.151 


0.206 


0.234 




3 P2 


0.200 








Ca 


'Pi 


0.216 


0.234 


0.211 


0.157 




3 Po 


0.138 










3 Pl 


0.139 


0.090 


0.144 


0.157 




3 P2 


0.140 








Zn 


l Pl 


0.426 


0.422 


0.403 


0.314 




3 Po 


0.294 










3 Pl 


0.296 


0.250 


0.316 


0.314 




3 P2 


0.300 








Sr 


'Pi 


0.198 


0.210 


0.193 


0.141 




3 Po 


0.130 










3 Pl 


0.132 


0.081 


0.135 


0.141 




3 P2 


0.136 








Cd 


'Pi 


0.398 


0.376 


0.346 


0.269 




3 Po 


0.274 










3 Pl 


0.279 


0.211 


0.272 


0.269 




3 P2 


0.290 









Table 2: The lowest S-^P excitation energies of various atoms. The experimental 
values (Erst column) [226] are compared with results calculated from Eq. (362) for 
the singlet and from Eq. (363) for the triplet (second column) and with ordinary 
Ascf values (third column). The optimized effective potential was used for u xc and 
the approximate OPM kernel (325) for / xc and G xc . The corresponding Kohn-Sham 
orbital-energy differences uj n are shown in the last column (All values in rydbergs). 
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Atom 


State 


^cxp 


^nPM-Li 1 1 1 A 
\.J r 1V1 -\- J\ Li D J\ 


^(Ascf) 


, 1Y1 

uj 


Be 


1 r> 
Pi 


0.388 


0.398 


0.331 


0.259 




Po 


0.200 










3 Pi 


0.200 


0.196 


0.181 


0.259 




3 P 2 


0.200 








A /Try 

IVlg 


1 P, 
"1 


u.oiy 


u.ozy 


n 9QQ 






3 Po 


0.199 










3 P1 


0.199 


0.196 


0.206 


0.234 




3 P2 


0.200 








Ca 


'Pi 


0.216 


0.236 


0.211 


0.157 




3 Po 


0.138 










3 Pl 


0.139 


0.129 


0.144 


0.157 




3 P2 


0.140 








Zn 


l Pl 


0.426 


0.417 


0.403 


0.314 




3 Po 


0.294 










3 Pl 


0.296 


0.280 


0.316 


0.314 




3 P2 


0.300 








Sr 


'Pi 


0.198 


0.211 


0.193 


0.141 




3 Po 


0.130 










3 Pl 


0.132 


0.117 


0.135 


0.141 




3 P2 


0.136 








Cd 


'Pi 


0.398 


0.370 


0.346 


0.269 




3 Po 


0.274 










3 Pl 


0.279 


0.239 


0.272 


0.269 




3 P2 


0.290 









Table 3: The lowest S-^P excitation energies of various atoms. The experimental 
values (Erst column) [226] are compared with results calculated from Eq. (362) for 
the singlet and from Eq. (363) for the triplet (second column) and with ordinary 
Ascf values (third column). The optimized effective potential was used for v xc 
and the ALDA for the xc kernels. The corresponding Kohn-Sham orbital-energy 
differences iv are shown in the last column (All values in rydbergs). 
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systems can be related to the electric polarizabilitics of the interacting species, 
and also to the correlation of fluctuating electric multipoles on the two systems. 
In the Present TDDFT context, a useful polarizability form for the second-order 
dispersion interaction was given by Zaremba and Kohn [231] who derived it directly 
from second-order perturbation theory: 

1 f ,Q f ,1 f ,1 / f ,1 i 1 1 



= -±- / d 3 n / d\ 2 / dM / a 



2wJ " -J ~ ''J ~ V " 2 In-rj 



x / duXa(ri,ir'i,iu)Xb(r2,r'2,iu)- (368) 
Jo 

Here Xa(r,r',w) and Xb( r , r ',^) are the exact density-density response functions 
(157) of each separate system in the absence of the other. Xa is defined by the 
linear density response ni a (r) exp(wi) of the electrons in system a to an externally 
applied electron potential energy perturbation V^ xt {v)e ut : 



ni a (r)= J dVxa^r'.i^y^fr') 



(369) 



and similarly for Xb- It is important to note that Xa includes the electron-electron 
interaction amongst the electrons of system a to all orders, and similarly for Xb- 
(Note also that, unlike Ref. [231], we have referred the space arguments of Xa and 
Xb in (368) to a common origin.) 

The expression (368) is more general than the familiar asymptotic i?~ 6 form. It 
applies to neutral quantal systems of any shape (not necessarily spherical) provided 
that R is still large enough that the electron densities do not overlap and that the 
inter-system Coulomb interaction can be treated in second order. We can recover the 
i?~ 6 form by assuming that R >> A, B where A and B are the spatial dimensions 
of the two systems. Then one can expand the Coulomb interactions in (368) in 
multipoles. The lowest nonvanishing term gives, with the "3" axis chosen along R, 

E^{R)~ ± t (1 ' 3 ^ ( a 1 ' 3 " 3j) I™ dua% { iu)^{iu), R » A, B 

i,j=l 

(370) 

where, for each system 

onM = j d\ j d 3 r' (n - Xi)(rJ - X,) X (r,r» (371) 

is the dipole polarizability tensor and X is the centre of electronic charge of the 
system. When the polarizabilities are isotropic so that ay = aSij, (370) reduces to 
the more familiar London form [229] 

Ei2) (R)~^J dua (a \iu)a^(iu). (372) 

Van Gisbergen, Snijders and Baerends [232] have evaluated a formula equivalent 
to (370) for diatomic and polyatomic molecules, using the ALDA to obtain the {a^ }. 
They find that, for the isotropic part of the vdW interaction, ALDA gives errors 
of similar size (but mostly opposite sign) to time-dependent Hartree Fock theory 
(except for the smallest atoms). This was achieved with much less computational 
effort than in the time-dependent Hartree Fock approach. The isotropic vdW coef- 
ficients, like the static and dynamic polarizabilities, were found to be somewhat too 
large. For the anisotropic part of the interaction, they found that ALDA compares 
favourably with both Hartree-Fock and Many-Body Perturbation Theory. Scincc 
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the ALDA contains unphysical orbital self-interaction, one can speculate that the 
use of self-interaction corrected (SIC) functionals might further improve the accu- 
racy of the method employed in Ref . [232] . These unphysical self- interactions cause 
orbitals to "see" an incorrectly large charge from the other electrons, causing or- 
bitals to be too spatially extended and hence too polarizable. This presumably has 
effects mainly at the isotropic level. Furthermore, SIC phenomena are known to be 
strongest for small systems with highly localized orbitals. Both of these consider- 
ations can be expected to cause difficulties in the very cases where van Gisbergen, 
Snijders and Baerends observed the least favourable ALDA results in comparison 
with other methods. This SIC explanation gains further support from the work of 
Pacheco and Ekardt [166] on alkali metal microclusters. Their static and dynamic 
SIC terms [233] were found to have significant effects on the polarizability and vdW 
interaction for small clusters and even, to a lesser degree, for quite large ones. 

Van Gisbergen et al. [232] commented that their numerical method could ac- 
commodate more sophisticated forms of TDDFT than simply the ALDA, and in 
particular, considering temporal and spatial nonlocality in the xc kernel, they felt 
that the latter might be the more important. 

Before leaving the discussion of vdW interactions in non-overlapping systems, 
we mention that the exact second-order dispersion formula (368) can be used [234] 
to derive a class of approximate vdW expressions for the groundstate energy as an 
explicit but highly nonlocal functional of the groundstate density. The idea is to 
make a direct local density approximation for the interacting susceptibilities Xa and 
Xb in (368). Extreme care is needed, however, to ensure one docs not violate the 
charge conservation condition 

|A X (r,r', W )= (373) 

or the reciprocity condition 

x(r, r', iu) = x(r', r, — iu) for real u. (374) 

An Ansatz satisfying these conditions and based on the simplest, pressure-free hy- 
drodynamic analysis of the uniform electron gas was given in [234]: 

X ^ o r(r,r',a;) = V r -V r - 

When this is substituted into (368) for each of Xa and Xb one obtains 

E {2) = -J^ I d s n [ d\ 2 ~ ^ , (376) 
327r 2 i V rf 2 (wx+wa) V ; 

where u>\ = topi = (47rn a (ri)/m) 1 / 2 is the local plasma frequency at an arbitrary 
point i"i inside system a, and similarly for uj 2 - Equation (376) constitutes a very 
nonlocal groundstate density functional, and it clearly provides a systematic basis 
for the much-used [229] simple notion of pairwise addition of i?~ 6 vdW contribu- 
tions. It is interesting that the integrand in (376) is proportional to the harmonic 
mean, + UJ2), of the two local plasma frequencies. The same formula (376) 

was very recently postulated [235] by Andersson, Langreth and Lundqvist on dif- 
ferent grounds. They obtained (376) by examining limiting cases and so modifying 
a somewhat similar formula previously postulated by Rapcewicz and Ashcroft [236] 
on the basis of diagrammatic arguments. The Rapcewicz-Ashcroft formula differs 
from (376) only in the replacement of oj\ + uj 2 by 2^/lJiuj2 on the denominator of 
(376). It was shown in [236] and [235] that these simple formulae give quite good 



n(r)<5(r — r ) 
lu 2 - uj 2 p Mr)) 



(375) 
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answers for the isotropic R~ 6 dispersion coefficient for various atomic pairs, pro- 
vided that one uses an appropriate cutoff in the low-density tails of the electron 
distributions. 

The derivation of (376) given in [234] promises to be extendable to more sophisti- 
cated local approximations for \a and Xb in (368), based perhaps on hydrodynamics 
with the inclusion of pressure (Thomas- Fermi hydrodynamics [112]) or of pressure 
plus density gradient (Thomas- Fermi- Weizsacker hydrodynamics [237] ) . With suit- 
able care to satisfy the constraints (373) and (374), one may thereby hope to obtain 
a more accurate extension of (376) involving gradients of the groundstate density 
and, possibly, having less dependence on spatial cutoffs. 

(ii) vdW interactions in closely juxtaposed or overlapping systems: The work of 
van Gisbergen et al. [232] and Pacheco and Ekardt [166], discussed in the previous 
section, shows that TDDFT, at least in the form of the ALDA, can represent the 
state of the art in evaluating van der Waals interactions in well-separated systems 
that are too large for methods such as the Configuration Interaction approach. 
What of more general cases where the electron clouds overlap or where no large 
separation exists? To study this for large systems, we seek a density functional 
approach, but first we need to appreciate the origin of the vdW force in terms of 
correlation physics. 

In essence, dispersion forces arise from the correlation between dynamic charge 
density fluctuations in two different systems or in distant parts of one system. The 
difficulty [228] in describing vdW forces in the static LDA or gradient approaches 
is therefore not surprising since in a highly inhomogeneous system (exemplified by, 
but not limited to, a pair of separated subsystems) these correlations may be quite 
different from those in the uniform or near-uniform electron gas upon which the 
LDA and the various gradient approximations are based. 

The previous section applied only to well-separated subsystems. The neces- 
sary correlations between distant fluctuations were generated by the application of 
second-order perturbation theory, and the TDLDA aspect of the calculation was 
not called upon to produce the vdW correlations directly. For overlapping systems 
(and for some closely juxtaposed systems), low-order perturbation theory in the 
Coulomb potential is not appropriate. The present section will outline an approach, 
currently under development, which does generate such long-ranged correlations in 
a natural fashion by the solution of a highly nonlocal real-space screening integral 
equation. Nevertheless, local density approximations are made wherever possible for 
the independent-electron susceptibility x.s an d the exchange-correlation kernel f xc , 
neither of which needs to be long-ranged in order to generate the basic long-ranged 
vdW correlations. 

The starting point for the proposed new approach is an exact formula [238] , [239] , 
based on the adiabatic connection formula and the zero-temperature fluctuation- 
dissipation theorem, relating the groundstate xc energy to the interacting suscepti- 
bility x ■ 

E XC = ~J\\ J d 3 r J d 3 r'j^-^ J™ du X (A, r, r', iu)j + n(r)5(r - r') 

(377) 

Here x(A, r, r', cu) is the interacting susceptibility defined as before but with a re- 
duced Coulomb interaction X/r acting between electrons. It was shown in Ref. [234] 
that the charge conservation condition (373) for \ implies xc hole normalization. 
Use of the independent-electron Kohn-Sham susceptibility Xs — x(A = 0,r,r',iu) 
rather than x(A, r, r', iu) in (377) yields the exact exchange energy. Subtraction 
of this exchange energy expression from the above xc energy yields the correlation 
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energy 

E c = -^J q dX J d3r J d3r 'y^\J du(x(A,r,r / ,i«)-Xa(r,r / ,i«)). (378) 

Equation (378) is required to produce the dispersion interactions under study. Pe- 
tcrsilka, Gossmann and Gross [69] have shown that \ an d Xs are related exactly by 
a Dyson-type equation involving the dynamic nonlocal xc kernel f xc as well as the 
Coulomb kernel (cf. Equation (152)): 

X(r,r',w) = x s (r,r',w)+ J d 3 x J d 3 x' Xs (r, x, w) ( j^~~7j + /xc(x,x',u)^ x(x',r',u) 

(379) 

Equations (377) and (379) are of course exact provided that f xc is exact, and 
so they contain inter alia the exact vdW interaction. Consider first a homogeneous 
electron gas. If f xc is arbitrarily set to zero, and (379) is Fourier-transformed with 
use of the convolution theorem, (379) is then seen to be the equation for the RPA 
response function \ m terms of the bare (dynamic Kohn-Sham-Lindhard) response 
Xs- Again with the assumption f xc = 0, but with the homogeneous assumption 
removed, (377) and (379) merely represent the inhomogeneous generalization of 
the well-known RPA groundstate correlation energy of the homogeneous electron 
gas. This case of zero /^already has some rather useful properties with respect to 
the vdW interaction. It has been shown in detail [239] that, when the correlation 
energy recipe (377), (379) with f xc = is applied to an arbitrary pair of widely- 
separated systems, the Zaremba-Kohn second-order vdW energy expression (368) is 
reproduced, with the following exception: the susceptibilities Xa and Xb are the ap- 
proximate RPA-interacting susceptibilities of each system, rather than including the 
exact interactions within each subsystem. Thus the full inhomogeneous RPA corre- 
lation energy already contains the essence of the vdW interaction, and will produce 
an R~ e dependence in the appropriate limit. An examination of the detailed proof 
in [239] further shows that the long-ranged vdW interaction achieved in the RPA 
does NOT arise because of any long-ranged behaviour of the independent-electron 
susceptibility Xs ( indeed Xs is not normally long-ranged). Rather, the long range 
of the vdW interaction comes from the long range of the Coulomb interaction in the 
screening equation (379). Thus a local density approximation for \s will not spoil 
the vdW properties, but may slightly alter the interacting susceptibilities Xa and Xb 
in the asymptotic form (368) . Furthermore the reintroduction of f xc within a local 
approximation can also be seen, from the working of Ref. [239], to maintain the 
form (368) in the separated limit, but the individual susceptibilities Xa and Xb will 
now involve f xc and hence will be closer to the required interacting susceptibilities. 

To summarize the previous paragraph: If we make short-ranged local-density 
or gradient approximations for Xs and f xc in the exact groundstate energy scheme 
(377), (379), we obtain an approximate and highly nonlocal prescription for the 
groundstate correlation energy, with the groundstate density n(r) as the only in- 
put. This scheme is expected to produce a rather good approximation to the long- 
ranged vdW dispersion interaction between widely separated subsystems, a result 
due principally to the retention, in full, of the nonlocal coulomb kernel in the real- 
space screening integral equation (379). 

What is now required is a sufficient set of constraints on the kernel f xc so that the 
short-ranged aspects of the groundstate correlation energy are also reproduced by 
(377), (379) at a level of approximation comparable, say, to the groundstate LDA or 
GGA. If this can be achieved, we will have a "seamless" scheme, equally reasonable 
for chemically bonded systems, metals etc., and also for fully or partly subdivided 
systems at all separations. This should allow investigation of the intermediate region 
of interaction where both short-ranged and long-ranged correlations are significant, 
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even for systems too large for traditionally accurate methods such as CI or IVfoller- 
Plesset perturbation theory. (Recall that wavefunctions are not needed for the 
present scheme, only groundstate densities, so that one may perform "real-space 
quantum chemistry" without basis-set problems). 

The details of this scheme are currently being worked out. Although it aims 
for a groundstate energy functional, it depends heavily on time-dependent density 
functional theory in the sense that the properties of the dynamic TDDFT xc ker- 
nel f xc (r,r',iu) for inhomogeneous systems are of the essence. Further details of 
some constraints to be obeyed by f xc are discussed in Ref. [239]. Some supporting 
evidence for the utility of a local approximation for \s m a highly inhomogeneous 
system are given in [113]. Current indications are that, for jellium slab situations 
where the inhomogeneity is only one-dimensional, the complete scheme (377), (379) 
(even with the exact Kohn-Sham \ s but with a local approximation for f xc ) can be 
computed on a single-processor 1- MFlop workstation in 10 2 hours or less. With 
Monte Carlo methods for the integrations in (377), and/or faster or parallel ma- 
chines, more involved geometries should be tractable. 

8 Applications beyond the perturbative regime: 
Atoms in strong femto-second laser pulses 

Owing to rapid experimental progress in the field of laser physics, ultra-short laser 
pulses of very high intensity have become available in recent years. The electric 
field produced in such pulses can reach or even exceed the strength of the static 
nuclear Coulomb field. If an atomic system is placed in the focus of such a laser 
pulse one observes a wealth of new phenomena [240] which cannot be explained by 
perturbation theory. In this case a non-perturbative treatment, i.e., the solution 
of the full TDKS equations (39) - (41) is mandatory. The total external potential 
seen by the electrons is given by 

v(r, t) = h E f(t)z sin(w t) (380) 

r 

where Z is the nuclear charge. The second term on the right-hand side of Eq. (380) 
is the potential due to the laser field in dipole approximation, written in the length 
form. Since the wavelength of currently used lasers is almost always very large 
compared to any characteristic length associated with an atomic system, the dipole 
approximation turns out to be very good in practice [241]. E denotes the peak 
electric field strength and f(t) characterizes the envelope function of the pulse which, 
in the calculations described below, is linearly ramped to its peak value over the 
first 10 cycles and then held constant. The field is assumed to be polarized along 
the z-direction. 

In the following, we compare the results of a TDKLI calculation using the ap- 
proximate potential (313) with an ALDA calculation using the potential (186), both 
for the exchange-only case [242, 243]. The numerical procedure [244] to solve the 
TDKS equations is similar to the one developed by Kulander [245, 246], who solved 
the time-dependent Schrddingcr equation for hydrogen and the time-dependent 
Hartree equation for helium in a laser pulse. The spin orbitals are expressed in 
cylindrical coordinates and, due to the linear polarization of the field, the spin 
as well as the angular part of the orbitals are preserved. Consequently, a fully 
three-dimensional treatment only requires a two-dimensional grid for the numerical 
integration. In the following, the time-dependent orbitals will always be character- 
ized by the indices indicating the initial state of the respective orbital; e.g., y2s(r, t) 
describes an electron which initially was in a 2s spin orbital: f2s( r , t = 0) = 02s ( r )- 
The integration of the single-particle equations is performed using a finite-difference 



65 



le-5 - 




Harmonic order 

Figure 4: Harmonic spectrum for He at A = 616nm and 7 = 3.5xl0 14 W/cm 2 . The 
squares represent experimental data taken from Ref. [248] normalized to the value 
of the 33rd harmonic of the calculated spectrum. The experiment was performed 
with a peak intensity of 1.4 x 10 14 W/cm 2 . 



representation of the kinetic energy operator. A Crank-Nicholson technique is em- 
ployed for the (unitary) time propagation of the orbitals. 

Once a numerical solution of the TDKS equations has been obtained, the result- 
ing time-dependent density is sufficient to calculate any desired observable of the 
system. Some quantities are easily calculated while others (such as ATI spectra) 
are harder to extract from the density. But, as demonstrated in section 2, all phys- 
ical observables can be calculated from the density, in principle. In the following 
we shall describe the calculation of two different quantities, namely the harmonic 
spectrum and the ionization yields. 

To obtain the harmonic spectrum, we calculate the induced dipole moment 

d(t) = J d 3 rzn(r,t) (381) 

which is then Fourier transformed over the last 5 cycles of the constant-intensity 
interval. The square of the resulting Fourier transform, |d(u)| 2 , has been shown 
[247] to be proportional to the experimentally observed harmonic distribution to 
within a very good approximation. Figure 4 shows the result of a simulation 
for the helium atom at a laser wavelength of A = 616 nm and peak intensity of 
I = 3.5 x 10 14 W/cm 2 . The calculation was made with the TDKLI scheme which, 
for two electrons in the x-only limit, reduces to the ordinary time-dependent Hartree 
method. One observes peaks in the energy-resolved photon spectrum at odd mul- 
tiples of the external laser frequency. From perturbation theory one would expect 
an exponential decrease of the peak intensities. Figure 4, however, shows a plateau 
of peak intensities up until roughly the 47th harmonic. This plateau is a typical 
nonlinear phenomenon. The squares in Fig. 4 indicate experimental results [248] 
obtained with the same laser frequency at an intensity of 1.4 x 10 14 W/cm 2 . Various 
calculations were performed with different peak intensities, but the best agreement 
with the experiment was achieved in the calculation for I = 3.5 x \Q 1A W/cm 2 shown 
in Fig. 4. The discrepancy between this intensity and the experimental intensity of 
1.4 x 10 14 W/cm 2 might be due to the uncertainty of the experimentally determined 
peak intensity which can be as high as a factor of two. 
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Figure 5: Harmonic distribution for He in a two-colour laser field. The two wave- 
lengths are 616 nm and 308 nm, and the intensity is 3.5 x 10 14 W/cm 2 for both of 
them. Crosses are the results for ip = and diamonds denote the values obtained 
with phase shift tp — 0.7ir. For comparison, the squares indicate the harmonic 
distribution for He in a one-colour field with A = 616 nm and I = 7 x 10 14 W/cm 2 . 



The harmonic generation of helium in a strong two-color laser field has also been 
studied [243, 108]. The two lasers with frequencies oj and 2w , respectively, are 
operated with the same peak intensity and a constant relative phase difference ip. 
This results in a total external potential of the form 

v(r, t) = h E f{t)z\sm(us t) + sin(2w t + (p)] (382) 

r 

where both fields are linearly polarized along the z-axis. 

Calculated harmonic distributions induced by a two-colour field with different 
relative phases are shown in Fig. 5. To avoid overcrowding, only the calculated peak 
intensities are plotted and connected with straight lines. The fundamental wave 
length is again 616 nm and the intensity is 3.5 x 10 14 W/cm 2 for both frequency 
components. We also show the one-colour spectrum for A = 616 nm calculated with 
the same total intensity as the two-colour field, i. e. I = 7 x 10 14 W/cm 2 . In the 
two-colour spectrum, harmonics at all higher multiples (including even multiples) 
of the fundamental frequency ujo occur due to nonlinear mixing processes of the two 
fields [249] . Most of the harmonics produced by the two-colour field in the plateau 
region are one to two orders of magnitude more intense than those obtained in 
the one-colour calculation although the total intensity of the external laser field is 
the same in all cases. Similar results have recently been found for hydrogen in a 
two-colour field [250]. One possible reason for this remarkable enhancement is that 
in a two-colour field one specific high-order harmonic can be generated by a large 
number of different mixing processes [249]. 

In order to simulate ionization, the grid contains an absorbing boundary to 
remove the flux of electrons leaving the nucleus. When some portion of the wave 
function propagates to the outer edges of the grid it is absorbed. We assume this 
flux corresponds to the ionized part of the wave function. Strictly speaking, such a 
criterion is meaningful only after long times when the respective contributions have 
propagated very far away from the nucleus. For the wave lengths considered here, a 
cylindrical grid of 20 x 60 a.u. was found to be sufficient. As time proceeds, more 
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and more electrons will be removed from the atom and, accordingly, the norm of 
the TDKS orbitals taken over the finite volume of the grid, 



decreases with time. 

Figure 6 compares the results of a TDKLI and an ALDA calculation [243] for 
Ne exposed to a laser field with wavelength A = 248 nm and intensity J = 3 x 
10 15 W/cm 2 . Figure 6 shows the norm (383) of those orbitals which were initially 
in the Ne 2s, 2po and 2pi states. The Is electrons have been frozen, i.e., only the 2s 
and 2p electrons are propagated by solving the TDKS equations, whereas the time 
evolution of the Is electrons is given by 



As expected, among the Ne 2s, 2p and 2p! orbitals, the 2s orbital is the least 
ionized one because it is initially more strongly bound (by roughly a factor of 2) 
than the 2p orbitals. A little surprising at first sight, the 2p and 2pi orbitals differ 
by about an order of magnitude in their degree of ionization (60% for the 2p orbital 
compared to only 4.75% for the 2pi orbital within TDKLI, and 56% for the 2po 
compared to 7.7% for the 2p! orbital within the ALDA). This difference has been 
observed before by Kulander [251, 252] for the case of xenon (in a single-active- 
electron calculation) . It is due to the fact that the 2po orbital is oriented along the 
polarization direction of the laser field, which makes it easier for the electrons to 
escape the nuclear attraction than for the case of the 2pi orbital, which is oriented 
perpendicularly to the field polarization. 

To explain the difference between the results obtained within the TDKLI and 
ALDA schemes shown in Fig. 6, we observe that the initial 2s and 2p , 2pi orbital 
energies in LDA differ quite considerably from those obtained with the KLI method: 
It takes 5 photons to ionize the 2p orbitals in KLI compared to only 3 photons in 
LDA. Similarly, it takes 11 photons to ionize the 2s orbital in KLI and only 9 in 
LDA. The difference between the curves in Fig. 6A and C is thus hardly surprising. 
On the other hand, it seems quite unexpected that the ALDA and TDKLI curves 
cross in Fig. 6B so that the ALDA curve comes to lie above the TDKLI curve. 
This behaviour can be attributed to the fact that the other orbitals are ionized 
much more strongly in ALDA than in TDKLI, so that their electron density near 
the nucleus (and therefore their screening of the nuclear charge) is decreased. This 
makes it more difficult for the 2p electrons to escape within the ALDA scheme. 

Figure 6 clearly shows the superiority of the TDKLI approach over the ALDA. 
The spurious self-interaction present in the ALDA causes the orbitals to be too 
weakly bound and hence the ALDA is not reliable in the calculation of ionization. 

The probabilities of finding neutral, singly, doubly, etc. ionized atoms at time 
t are readily expressed in terms of the norms (383). For instance, in the case of 
helium, one has 
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and the probabilities for neutral, singly and doubly charged helium are 



P°(t) = N ls (t) 2 
P +1 (t) = 2JVi„(t)(l- N ls (t)) 
P +2 (t) = (l-N ls (t)) 2 . 



(386) 
(387) 
(388) 
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Figure 6: Time evolution of the norm of the Ne 2s orbital (A), the Ne 2po orbital 
(B) and the Ne 2pi orbital (C), calculated in the x-only TDKLI and ALDA schemes. 
Laser parameters: A = 248 nm, J = 3 x 10 15 W/cm 2 , linear ramp over the first 10 
cycles. One optical cycle corresponds to 0.82 femtoseconds. 
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Figure 7: Population of the differently charged states of Ne. Laser parameters as 
in Fig. 6 (X — 248 nm, I = 3 x 10 15 W/cm 2 , linear ramp over the first 10 cycles). 

For many-electron atoms similar combinatorical considerations [244] are performed 
to determine the probabilities for the various charged ions. Figure 7 shows the 
probabilities of finding neutral, singly, doubly and triply charged Ne as calculated 
from the norms of Fig. 6. 

These probabilities as a function of time cannot be compared directly with 
experiment. This is because the laser focus, in addition to the temporal pulse 
shape, has a spatial intensity profile due to which not all atoms in the laser focus 
experience the same intensity. Hence a realistic calculation of ion yields requires 
many runs at various peak intensities. Work along these lines remains an important 
field for the future. In this way one might be able to understand the structures in 
the strong field ionization spectra of He [253] which have been the subject of heated 
discussions in recent years. 
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